QUASI-TRANSIENT  CALCULATION  OF  SURFACE  TEMPERATURES  ON  A  REUSABLE 


BOOSTER  SYSTEM  WITH  HIGH  ANGLES  OF  ATTACK 


Thesis 

Submitted  to 

The  School  of  Engineering  of  the 
UNIVERSITY  OF  DAYTON 


In  Partial  Fulfillment  of  the  Requirements  for 
The  Degree  of 

Master  of  Science  in  Aerospace  Engineering 


By 

Seth  Henderson  Morris 
Dayton,  OH 
December,  2011 


UNIVERSITY  of 


DAYTON 


Form  Approved 
OMB  No.  0704-0188 


Report  Documentation  Page 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

DEC  2011 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2011  to  00-00-2011 

4.  TITLE  AND  SUBTITLE 

Quasi-Transient  Calculation  of  Surface  Temperatures  on  a  Reusable 
Booster  System  with  High  Angles  of  Attack 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Dayton, 300  College  Park, Dayton, OH, 45469 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS (ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

The  calculation  of  a  recovery  temperature  based  heat  transfer  coefficient  proves  to  be  sufficiently 
independent  of  wall  temperature  to  use  in  a  three  dimensional,  transient  temperature  model  of  a  thermal 
protection  system  of  a  reusable  booster  concept.  After  a  derivation  of  recovery  temperature  from  the  1st 
law  of  thermodynamics,  the  weak  dependence  of  the  recovery  temperature  based  heat  transfer  coefficient 
is  investigated  by  72  Computational  Fluid  Dynamics  (CFD)  models  at  angles  of  attack  ranging  from  0?  to 
90?  over  a  range  of  Mach  numbers,  from  Mach  2  to  5,  and  a  variety  of  thermal  boundary  conditions  at  the 
wall,  from  isothermal  to  a  conductive  wall.  Then,  the  heat  transfer  coefficient  is  calculated  at  many  steady 
state  CFD  solutions  for  a  reusable  booster  system  concept  on  a  given  trajectory  and  applied  to  a  transient 
Finite  Element  Analysis  (FEA)  model  of  a  thermal  protection  system.  Results  are  presented  graphically. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

110 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


QUASI-TRANSIENT  CALCULATION  OF  SURFACE  TEMPERATURES  ON  A  REUSABLE 
BOOSTER  SYSTEM  WITH  HIGH  ANGLES  OF  ATTACK 


Name:  Morris,  Seth  Henderson 


APPROVED  BY: 


Timothy  John  Fry,  Ph.D. 

Advisory  Committee  Chairman 
Group  Leader 

Experimental  and  Applied  Mechanics  Group, 
University  of  Dayton  Research  Institute 


Jose  A.  Camberos,  Ph.D.,  P.E. 

Committee  Member 
Professor 

Department  of  Mechanical  and  Aerospace 
Engineering 


John  Doty,  Ph.D. 

Committee  Member 
Professor 

Department  of  Engineering  Management 
and  Systems 


John  G.  Weber,  Ph.D. 
Associate  Dean 
School  of  Engineering 


Tony  E.  Saliba,  Ph.D. 

Dean,  School  of  Engineering 
&  Wilke  Distinguished  Professor 


ABSTRACT 


QUASI-TRANSIENT  CALCULATION  OF  SURFACE  TEMPERATURES  ON  A  REUSABLE 
BOOSTER  SYSTEM  WITH  HIGH  ANGLES  OF  ATTACK 

Name:  Morris,  Seth  Henderson 
University  of  Dayton 

Advisor:  Dr.  Timothy  J.  Fry 

The  calculation  of  a  recovery  temperature  based  heat  transfer  coefficient  proves  to 
be  sufficiently  independent  of  wall  temperature  to  use  in  a  three  dimensional,  transient 
temperature  model  of  a  thermal  protection  system  of  a  reusable  booster  concept.  After  a 
derivation  of  recovery  temperature  from  the  1st  law  of  thermodynamics,  the  weak 
dependence  of  the  recovery  temperature  based  heat  transfer  coefficient  is  investigated  by 
72  Computational  Fluid  Dynamics  (CFD)  models  at  angles  of  attack  ranging  from  0°  to  90° 
over  a  range  of  Mach  numbers,  from  Mach  2  to  5,  and  a  variety  of  thermal  boundary 
conditions  at  the  wall,  from  isothermal  to  a  conductive  wall.  Then,  the  heat  transfer 
coefficient  is  calculated  at  many  steady  state  CFD  solutions  for  a  reusable  booster  system 
concept  on  a  given  trajectory  and  applied  to  a  transient  Finite  Element  Analysis  (FEA) 
model  of  a  thermal  protection  system.  Results  are  presented  graphically. 
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CHAPTER  1 


INTRODUCTION  AND  PROBLEM  DESCRIPTION 


The  calculation  of  external  surface  temperatures  due  to  aerodynamic  heating  often 
presents  technical  difficulties  for  thermal  protection  system  (TPS)  engineers.  The  transient 
surface  temperatures  of  a  high  speed  vehicle  are  difficult  to  calculate  due  to  numerous 
factors.  Examples  of  these  difficulties  are  transient  flow-field  effects,  material  properties  of 
the  TPS,  and  accuracy  of  the  solver  or  the  solution  process.  The  complexities  involved  in  the 
calculation  of  the  transient  temperatures  present  a  need  for  research  of  improved  methods 
of  analysis.  In  the  past,  TPS  temperature  calculation  methods  relied  on  highly  simplified 
models.  With  modern  computing  power,  these  models  have  the  potential  to  be  more 
realistic;  thus,  the  research  for  these  alternate  methods  of  calculation  tends  towards  higher 
fidelity.  This  thesis  investigates  a  process  for  calculating  transient  surface  temperatures  on 
a  hypersonic  vehicle  that,  at  some  points  in  a  given  flight  trajectory,  flies  at  high  angles  of 
attack.  The  improved  analysis  method  described  in  this  thesis  utilizes  Navier-Stokes 
Computational  Fluid  Dynamics  (CFD)  while  decoupling  the  aerothermal  heating  from  the 
transient  temperature  profile  of  the  TPS. 

Past  methods  of  calculation  of  temperature  profiles  used  in  sizing  the  TPS  include 
recovery  temperature  calculations  for  peak  heating  at  the  leading  edges,  lower  fidelity, 
engineering  level  codes  such  as  Aeroheating  and  Thermal  Analysis  Code  (ATAC)  and 
Hypersonic  Engineering  Aerothermodynamic  Trajectory  Tool  Kit  (HEAT-TK),  and  steady 
state  CFD.  The  ultimate  model  of  heat  transfer  is  a  transient,  conjugate,  coupled  heat 
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transfer  model  where  the  internal  structure,  heat  sources  and  sinks,  TPS,  rocket  engines, 
and  plumes  are  all  modeled  along  with  a  transient  external  aerothermal  calculation. 
However,  this  level  of  completeness  and  integration  is  not  yet  possible  in  a  rapid  design 
environment  as  the  time-steps  required  for  external,  turbulent,  separating,  hypersonic  flow 
are  many  orders  of  magnitude  smaller  than  the  length  of  the  flight  trajectory.  For  this 
reason  we  seek  a  method  that  decouples  the  external  aerothermal  calculation  from  the 
transient  TPS  temperature  calculation. 

The  decoupled  method  performed  here  uses  steady  state  CFD  solutions  at  specific 
points  in  the  trajectory  to  generate  a  field  of  heat  transfer  coefficients  based  on  the  recovery 
temperature: 


In  Equation  (1-1)  q  is  the  heat  flux  per  unit  area,  Tw  is  the  wall  temperature,  and  Tr  is  the 
recovery  temperature.  Numerous  heat  transfer  coefficient  fields  are  calculated  over  the 
entire  vehicle  surface  based  on  far  field  pressure,  far  field  temperature,  free  stream  Mach 
number,  and  vehicle  angle  of  attack.  This  collection  of  fields  covers  the  full  operating  range 
of  the  vehicle  for  the  trajectories  to  be  analyzed.  Next,  a  particular  field  from  the  collection 
is  selected  based  on  the  flight  conditions  of  the  trajectory  time-step  under  consideration. 
This  particular  filed  is  held  constant  throughout  the  time-step  and  a  different  field  is 
selected  for  the  subsequent  time-step.  The  specific  heat  transfer  coefficient  field  serves  as 
the  convective  boundary  condition  on  a  TPS  transient  thermal  model.  This  process  is 
repeated  for  each  time-step  in  the  entire  trajectory  resulting  in  a  complete  calculation  of  the 
transient  TPS  temperature  profile. 

In  order  to  employ  this  method,  the  variation  associated  with  the  recovery 
temperature  based  heat  transfer  coefficient  with  regard  to  wall  temperature  must  be 
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investigated.  It  is  well  known  that  a  heat  transfer  coefficient  calculated  with  respect  to 
either  free  stream  or  boundary  layer  edge  temperature  is  highly  dependent  on  wall 
temperature  variations.  However,  a  HTC  based  on  recovery  temperature  does  not  exhibit 
this  tendency.  For  this  reason,  the  recovery  temperature  based  heat  transfer  coefficient  is 
significantly  more  useful  for  solving  the  TPS  temperature  problem  in  a  decoupled  manner. 
There  is  a  large  amount  of  work  performed  in  quantifying  the  wall  temperature  dependency 
of  the  heat  transfer  coefficients  based  on  the  thermal  boundary  conditions  of  the  wall  used 
in  the  solution  process. 

Quantifying  the  variation  of  the  recovery  temperature  based  heat  transfer  coefficient 
with  regards  to  wall  temperature  follows  this  process: 

1.  Generate  CFD  data  with  various  thermal  boundary  conditions  at  the  wall. 

2.  Calculate  the  recovery  temperature  at  each  CFD  cell  along  the  wall  for  each 
boundary  condition. 

3.  Calculate  the  heat  transfer  coefficient  based  on  recovery  temperature,  Equation 
(1-1).  This  is  performed  for  each  wall  boundary  condition. 

4.  Calculate  the  mean  heat  transfer  coefficient  of  all  the  wall  boundary  conditions. 

5.  Calculate  the  percent  difference  from  the  mean. 

6.  Investigate  the  results 


Figure  1-1  -  Geometry  of  NACA  23012  airfoil. 
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In  order  to  rapidly  quantify  the  wall  temperature  dependency  of  the  heat  transfer 
coefficient  based  on  recovery  temperature,  a  simple  geometry  is  used.  The  geometry  chosen 
is  a  NACA  23012  airfoil  with  a  2.5m  chord  length,  as  depicted  in  Figure  1-1.  The  CFD  data  is 
generated  at  various  flight  conditions,  such  as  different  angles  of  attack  and  Mach  numbers. 


Figure  1-2  -  Geometry  of  a  reusable  booster  system  concept. 

Once  the  variability  in  the  recovery  temperature  based  heat  transfer  coefficient  is 
quantified,  this  method  is  applied  to  a  three-dimensional  geometry  flying  a  given  trajectory. 
The  geometry  used  is  a  reusable  booster  system  concept  supplied  by  the  United  States  Air 
Force  Research  Labs,  shown  in  Figure  1-2.  The  trajectory,  internally  supplied,  shows  three 
stages:  ascent,  shown  in  blue,  rocket-back,  shown  in  red,  and  entry,  shown  in  green.  The 
vehicle  in  this  trajectory  flies  at  a  maximum  speed  of  just  over  Mach  6.3,  a  maximum  angle 
of  attack  of  approximately  151°,  and  a  maximum  altitude  of  just  over  300,000  feet. 
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Figure  1-3  -  The  trajectory  used  for  a  reusable  booster  system  concept. 
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CHAPTER  2 


THE  DERIVATION  OF  RECOVERY  TEMPERATURE  BASED  HEAT  TRANSFER  COEFFICIENT 
FROM  THE  FIRST  LAW  OF  THERMODYNAMICS 


In  this  chapter,  the  concepts  are  developed  in  order  to  calculate  temperature  rise  in 
a  fluid  due  to  a  rapid  change  in  velocity.  Working  in  a  vehicle-centered  coordinate  system, 
the  fluid  experiences  this  rapid  deceleration  when  it  flows  from  the  far-field  domain  to  the 
wall  of  the  vehicle,  where  the  fluid  reaches  zero  velocity.  Finding  an  expression  for  this  rise 
in  temperature  is  best  done  using  thermodynamic  processes,  starting  from  the  first  law  of 
thermodynamics. 

We  are  concerned  with  high  speed,  ideal  gas  flows.  Because  the  kinetic  energy  is  so 
dominant  in  this  type  of  flow,  the  gravity  vector  can  be  ignored;  thus,  potential  energy  can 
be  ignored. 


The  First  Law  of  Thermodynamics 

The  first  law  of  thermodynamics:  “during  any  cycle  a  system  (control  mass) 
undergoes,  the  cyclic  integral  of  the  heat  is  proportional  to  the  cyclic  integral  of  the  work" 
(Sonntag,  Borgnakke  and  Van  Wylen  116).  Mathematically,  then,  this  is  represented  as  the 
following  equation: 


(2-1) 
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If  we  take  Equation  (2-1)  in  the  context  of  a  cycle  that  goes  from  state  1  to  state  2,  and  use 
energy,  E,  we  can  analyze  Equation  (2-1).  Graphically,  we  have  illustrated  this  in  Figure  2-1. 


- ►  V 

Figure  2-1  -  A  fully  reversible  thermodynamic  process. 


The  mathematical  representation  of  Figure  2-1,  for  processes  A  and  B  is  as  follows: 

2  12  1 

1  SQA  +  j  SQb  =  I  SWA  +  I  SWB  (2-2) 

12  12 

Processes  B  and  C  are  as  follows: 

2  12  1 

jsQc  +  j  sQb  =  j  SWC  +  j  SWB  (2-3) 

12  12 

Subtracting  Equation  (2-3)  from  Equation  (2-2),  we  have: 

2  2  2  2 

j  SQa~  j  SQC  =  j  SWA-j  SM/C 
1111 


Or, 


2  2 

j  ( SQ  -  8W)a  =  J  (8Q  —  8W)C 
1  1 


(2-4) 
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Because  A  and  C  represent  arbitrary  processes  between  states  1  and  2,  the  quantity 
SQ  —  SW  is  the  same  for  all  processes  between  states  1  and  2.  Therefore,  SQ  —  SW  depends 
only  on  the  initial  and  final  states  and  not  on  the  path  followed  between  the  two  states 
(Sonntag,  Borgnakke  and  Van  Wylen  118).  Hence,  we  deduce  that  the  process  is  really  a 
“point  function,"  or  it  can  be  represented  as  a  differential  property  of  the  mass.  This 
differential  property  is  called  energy: 

dE  —  SQ  —  SW  (2-5) 

Integrating  Equation  (2-5)  from  point  1  to  2: 

E2-E1=1Q2-1W2  (2-6) 

Where  Et  and  E2  are  the  initial  and  final  values  of  the  energy  in  the  control  mass,  Q  is  the 
heat  transferred  to  the  control  mass  during  the  process  from  state  1  and  2,  and  W  is  the 
work  done  by  the  control  mass  during  the  process.  The  energy  term  may  be  broken  down 
into  three  basic  components: 

E  =  U  +  KE  +  PE  (2-7) 

Where  U  is  the  internal  energy,  KE  is  the  kinetic  energy,  and  PE  is  the  potential  energy. 
Equation  (2-7)  may  also  be  written  in  differential  form: 

dE  =  dU  +  d{KE )  +  d(PE)  (2-8) 

Combining  Equations  (2-8)  and  (2-5),  we  have: 

dU  +  d{KE )  +  d(PE)  =  SQ  -  SW  (2-9) 

In  order  to  find  the  kinetic  energy,  we  assume  no  heat  transfer  or  change  in  internal 
energy.  Recalling  that  W  =  F  ■  x,  Equation  (2-9)  simplifies  to: 

SW  =  -F  ■  dx  =  -d{KE ) 
dv  dx 

F  =  ma  —  m—  =  m  —  (V  ■  v)  =  rnv(V  ■  v) 
dt  dt 

Combining  the  above  two  equations,  we  have: 
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d(KE )  =  F  ■  dx  =  mv  ■  dv 


Integrating  the  above  equation,  we  get: 


d(KE)  =  I  mv  •  dv 


KE  =  -mv2 
2 


(2-10) 


We  will  define  the  kinetic  energy  per  unit  mass  as 


ke  =  -v2 
2 


(2-11) 


Equation  (2-11)  will  become  useful  later,  in  our  discussion  of  total  enthalpy. 

In  order  to  express  internal  energy,  we  substitute  the  differential  expression  for  kinetic 
energy  into  Equation  (2-8): 


dE  =  dU  +  mv  ■  dv 


(2-12) 


Integrating  the  above  equation, 


E2  —  Ei  —  U2  —  Ui  + 


mv,  mv! 


(2-13) 


If  we  take  the  differential  expression  for  kinetic  energy  and  substitute  it  into  Equation  (2-9), 
we  have: 


d(mv2) 

dU  +  „  =  SQ  -  SW 

2 


(2-14) 


Substituting  Equation  (2-13)  into  (2-6), 


(U2  -  Ut)  + 


mv2  mv{ 
~~2  2~ 


—  l(?2— 1^2 


(2-15) 


Equation  (2-15)  is  the  form  of  the  energy  equation  that  will  be  used  here. 
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Enthalpy,  Internal  Energy,  and  Specific  Heat 

Enthalpy  is  defined  by  considering  a  system  that  has  no  change  in  kinetic  energy. 
When  considering  this  type  of  system,  Equation  (2-15)  reduces  to: 

(U2-U1)=1Q2-1W2 

And  the  work  is  simply 

2 

J  PdV  =  tW2 
1 

Because  the  pressure  is  held  constant  in  this  example,  it  falls  outside  the  integral,  and 

P(V2  -  VO  =  xw2 

And 


(U2+P2V2)-(U1  +  P1V1)=1Q2 

Thus,  enthalpy  can  be  defined  by  the  following  equation  (Sonntag,  Borgnakke  and  Van 
Wylen  130): 

H  =  U  +  PV 

And  on  a  per  unit  mass  basis,  the  specific  enthalpy  h  is: 

h  =  u  +  Pv  (2-16) 


If  we  do  not  neglect  kinetic  energy,  which  is  necessary  high-speed  flows,  then  Equation  (2- 
15)  becomes: 

,  .  /mvf  mv}\ 

(y2-t/i)  +  hr  — A  =  1Q2-1W2 


Which,  using  the  same  analysis  as  above,  we  have: 


(  mv2\  (  ravA 

(  U2  +  P2V2  +  -  (A  +  PXVX  +  -y-J  =  tQ2 


The  above  equation  can  be  used  to  define  stagnation  enthalpy,  which  accounts  for  the 
kinetic  energy  of  the  fluid: 
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1  - 

HS  =  U  +PV  +  -mv2 
2 

The  above  equation  in  a  per-unit  mass  basis: 

1  _ 

h,  =  u  +  Pv  +  -  v 
2 


Or  simply, 


K  =  h  +  -v2 


(2-17) 


Equation  (2-17)  is  simply  the  sum  of  the  specific  enthalpy,  (2-16)  and  the  specific  kinetic 
energy,  (2-11). 

The  constant  pressure  specific  heat,  cp,  is  frequently  used  in  thermodynamic 
calculations.  If  we  start  with  Equation  (2-14),  and  disregarding  the  term  for  kinetic  energy, 
we  have: 


SQ  =  SW  +  dU 
SQ  =  dU  +  PdV 

Integrating  the  above  equation: 

2  2  2 

j  dQ  =  j  dU  +  j  PdV 
111 

P  is  constant  in  this  case,  so  we  have: 

Q2  ~  Qi  =  u2  ~  ux  +  P(v2  -  vi) 

The  resulting  P(V2  —  Vi)  term  can  be  associated  with  the  internal  energy  terms, 
similarly  to  enthalpy,  which  leads  to  the  conclusion  that  the  heat  transfer  can  be  expressed 
in  terms  of  the  enthalpy  change.  This  results  in: 


_  1  /SQ\  _  1  /dH\  _  /dh\ 
Cp  ~m\Sf)p  ~  rn\df)p  ~  Vd77P 

At  constant  volume,  the  work  term  PdV  is  zero.  This  results  in: 


(2-18) 
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(2-19) 


_  1  (SQ\  _  1  (dU\  _  (du\ 

Cv~m  \57/v  “  m  \df)v  ~  V  df)v 

Although  the  derivations  of  specific  heat  were  found  by  a  specific  process  that  neglected 
certain  thermodynamic  aspects  such  as  kinetic  energy,  the  resulting  definition  of  specific 
heat  is  path  independent: 

Note  that  in  each  of  these  special  cases,  the  resulting  expression,  [Equations  (2-18) 
and  (2-19)],  contains  only  thermodynamic  properties,  from  which  we  conclude  that 
the  constant  volume  and  constant  pressure  specific  heats  must  themselves  be 
thermodynamic  properties.  This  means  that,  although  we  began  this  discussion  by 
considering  the  amount  of  heat  transfer  required  to  cause  a  unit  temperature 
change  and  then  proceeded  through  a  very  specific  development  leading  to 
[Equations  (2-18)  and  (2-19)],  the  result  ultimately  expresses  a  relation  among  a 
set  of  thermodynamic  properties  and  therefore  constitutes  a  definition  that  is 
independent  of  the  particular  process  leading  to  it  (in  the  same  sense  that  the 
definition  of  enthalpy  is  independent  of  the  process  used  to  illustrate  one  situation 
in  which  the  property  is  useful  in  a  thermodynamic  analysis).  (Sonntag,  Borgnakke 
and  Van  Wylen  134) 

Recalling  the  ideal  gas  equation, 


Pv  =  RT  (2-20) 

And  the  definition  for  specific  enthalpy,  Equation  (2-16): 

h  =  u  +  Pv 

Combining  Equations  (2-20)  and  (2-16),  we  have: 

h-u  +  RT 

If  the  internal  energy  is  a  function  of  temperature  only,  u  =  f(T),  then  the  enthalpy  is  a 
function  of  temperature  only,  h  =  f(T ).  Thus,  the  partial  derivatives  in  Equations  (2-18) 
and  (2-19)  are  actually  ordinary  derivatives  for  ideal  gasses: 


c 


v 


(2-21) 


And 
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(2-22) 


cv  — 


v 


Because  h-u  +  jpv-u-\-  RT,  then 


dh  dv  dT 
dJ  =  dT  +  RdT 


And 


dh  =  dv  +  RdT  (2-23) 

Solving  for  dh  in  Equation  (2-21)  and  du  in  (2-22),  and  substituting  these  terms  into 
Equation  (2-23): 


cpdT  =  cvdT  +  RdT 


Collecting  terms,  and  integrating  with  respect  to  temperature, 

rT2 

0=1  (Cv  +  R  —  Cp)  dT 
0  =  (cv  +  R  —  cp)(T2  —  rx) 


Thus, 


Cy  T  R  cv 


(2-24) 


We  can  also  define  the  ratio  of  specific  heats,  y  =  —  in  an  alternate  fashion,  using  Equation 


(2-24): 


R 

y  =  i  +  — 


Or 


1  + 


=  Y 


cp  -  R  +  R  =  y(cp  -  R) 


cp  =  ycv  -  Ry 
ycp  —  cp  =  Ry 
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cp(y-  i)  =  Ry 


(2-25) 


Equation  (2-25)  will  become  useful  later. 

Note  that  in  our  discussion  leading  to  Equations  (2-21)  and  (2-22),  we  assume  an 
ideal  gas,  which,  at  the  molecular  level,  assumes  a  limited  number  of  degrees  of  freedom. 
For  example,  a  diatomic  molecule,  such  as  air,  can  be  modeled  “as  a  rigid  dumbbell,  with 
[number  of  degrees  of  freedom]  of  5“  (Liepmann  and  Roshko  36).  The  specific  heat  of  these 
diatomic  gasses  follow  their  ideal  gas  assumption  over  a  wide  range  of  temperatures,  but  at 
high  temperatures,  cp  falls  away  from  this  assumption  “because  the  atoms  in  a  diatomic 
molecule  are  not  rigidly  bound  but  can  vibrate,"  which  classically  adds  two  degrees  of 
freedom  (Liepmann  and  Roshko  36). 

Real  gasses  are  hard  to  model  exactly,  but  many  experiments  have  drawn  empirical 
relations  which  can  closely  approximate  a  real  gas's  behavior.  According  to  Liepmann  and 
Roshko,  “the  specific  heat  for  a  diatomic  gas  in  the  range  of  temperatures  between 
dissociation  and  the  region  of  van  der  Waals'  effects  is  closely  approximated  by: 


where  9V  denotes  a  characteristic  temperature  for  the  vibrational  energy"  (Liepmann  and 
Roshko  37). 

In  the  CFD  analyses  in  the  following  chapters,  the  peak  temperatures  in  the  flow- 
field  never  reach  the  temperatures  required  for  dissociation  of  the  molecules  that  comprise 
air.  Hirschel  provides  a  graph  that  relates  mole  percentage  of  different  materials  in  air 
versus  the  temperature,  Figure  2-2.  These  materials  are  diatomic  nitrogen,  N2 ;  diatomic 
oxygen,  02 ;  carbon  dioxide,  C02 ;  carbon  monoxide,  CO;  carbon,  C;  nitric  oxide,  NO; 
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monatomic  nitrogen,  N;  monatomic  oxygen,  0 ;  Argon,  Ar ;  electrons,  e ;  ionized  monatomic 
nitrogen,  N+;  ionized  monatomic  oxygen,  0+ ;  ionized  argon,  Ar+;  and  ionized  nitric  oxide 
NO+. 


T,  1000  "K 

Figure  2-2  -  Equilibrium  composition  of  air  at  p  =  10-2  atm  as  a  function  of  the 
temperature  (Hirschel  22). 

We  can  see  that  in  Figure  2-2,  the  real  gas  effects,  which  come  with  the  ionization 
and  dissociation  of  the  particular  molecules  that  comprise  air,  do  not  begin  to  occur  until 
the  temperature  reaches  ~3,500  Kf  where  diatomic  nitrogen  dissociates  to  two  monatomic 
nitrogens.  The  mole  percentages  of  monatomic  nitrogen  are  quite  small,  <  1%,  up  to  ~4,500 
K  The  approximation  of  an  ideal  gas  for  this  thesis  is  therefore  valid  as  the  temperatures  of 
the  flow-field  never  go  much  above  2,000  K 


15 


The  Speed  of  Sound 


P 

T 


P 

a 


Figure  2-3  -  Flow  through  a  disturbance. 


P  +  dP 
T  +  dT 


p  +  dp 
CL  +  dCL 


Consider  Figure  2-3,  where  P  is  the  pressure,  T  is  the  temperature,  p  is  the  density, 
and  a  is  the  speed  of  sound.  From  the  continuity  equation,  we  have: 


PlVl  =  P2V2 

Applying  Equation  (2-26)  to  Figure  2-3  results  in: 

pa  =  (p  +  dp)  (a  +  da) 

Expanding,  and  neglecting  terms  of  second-order  smallness, 

o 

pa  =  pa  +  pda  +  adp  +  dpda 
Reducing,  and  solving  for  the  speed  of  sound: 

da 


(2-26) 


(2-27) 


From  the  1-dimensional  momentum  equation: 

Pi  +  PiVi  =P2+  P2v| 

Applying  Equation  (2-28)  to  Figure  2-3,  we  have: 

P  +  pa2  =  (P  +  dP)  +  (p  +  dp)  (a  +  da)2 

Again,  expanding,  and  using  second-order  smallness  for  the  product  of  differentials, 


(2-28) 


P  +  pa2  =  (P  +  dP)  +  (p  +  dp)  a2  +  2 ada  +  (da)2 
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0 


P  +  pa 2  =  P  +  dP  +  pa 2  +  2  pada  4-  a2dp  +  2  adpda 
P  +  pa2  =  P  +  dp  +  pa2  +  a2dp  +  2  pada 
0  =  dp  4-  a2dp  +  2 pada 
2  pada  =  —dp  —  a2dp 


da  — 


—dp  —  a2dp 
2  pa 


Solving  Equation  (2-27)  for  da  and  substituting  the  resulting  term  into  (2-29), 


/-dp  -  a2dp\  _p  dp  +  a 
\  2  padp  J  2  ap 


(2-29) 


Solving  for  a, 


a 


2 


1 

2 


(dP  ,  2 

l TP  +  a 


) 


2  a2  —  a 2 


dp 

dp 


a2  =-f-  (2-30) 

dp 

Because  flow  through  a  sound  wave  is  isentropic,  Equation  (2-30)  is  the  rate  of 
change  of  pressure  with  respect  to  density,  and  the  ordinary  derivative,  which  is  an 
isentropic  change,  can  be  expressed  as  a  partial  derivative.  Hence, 


a  = 


(2-31) 


If  we  assume  an  ideal  gas,  the  isentropic  relation  given  by  the  following  Equation  holds 
(Anderson  558): 


Pi 

P2 


(2-32) 


From  Equation  (2-32),  we  get: 
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(2-33) 


Where  c  is  a  constant.  Solving  Equation  (2-33)  for  P  and  differentiating  with  respect  to  p, 
we  obtain: 


TP)S  =  w 


(2-34) 


Substituting  the  expression  for  the  constant,  Equation  (2-33),  into  (2-34): 


dP\  P 


rP>rprp 


dP  _  yP 
dp  p 


(2-35) 


Substituting  Equation  (2-35)  into  Equation  (2-31): 


If  we  substitute  the  ideal  gas  law,  -  =  RT  into  Equation  (2-36),  we  have: 


(2-36) 


=  JyRT 


(2-37) 


a2  =  yRT 


(2-38) 


Recovery  Temperature 

Using  the  relationship  expressed  in  Equation  (2-38)  and  taking  a  ratio  of  the  speed 

of  sound  at  the  stagnation  point  and  the  speed  of  sound  at  an  arbitrary  point,  we  see  that: 

qf  _  yRTs 
a2  yRT 


a£  =  Ts_ 
a2  ~  T 


(2-39) 
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Recall  the  derivative  definitions  of  the  specific  heats,  Equations  (2-21)  and  (2-22).  Solving 
for  the  differential  of  enthalpy  and  internal  energy,  we  have: 

dh  =  cpdT  (2-40) 

du  =  cvdT  (2-41) 

We  can  integrate  Equation  (2-40): 

h  =  cpT  (2-42) 

And  recalling  the  definition  for  stagnation  enthalpy,  we  can  combine  Equations  (2-17)  and 
(2-42): 

1 

hs  =  cpT  +  -v2 

Because  hs  =  cpTs,  we  have: 

1  _ 

cpTs  =  cpT  +  -v2 
Or 


T,  =  T  + 


2  c„ 


(2-43) 


If  we  take  the  ratio  of  the  stagnation  temperature  Ts  and  the  temperature  at  some 
arbitrary  point  T  in  Equation  (2-43),  we  arrive  to  the  following  expression: 

Ts  v2 

—  =  1  + - 

T  2  cpT 

From  our  relationship  found  in  Equation  (2-25),  we  can  revise  the  above  expression: 

v2(y- 1) 

T  ~  2 yRT 

Also,  recall  Equation  (2-38),  yRT  =  a2.  The  above  Equation  then  becomes: 

Ts  v2(y-  1) 

T  2a2 


Note  the  term  — ,  which  is  equal  to  the  square  of  the  Mach  number.  Hence, 
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(y-D 

2 


M2 


At  last,  we  have  an  equation  for  stagnation  temperature: 


TS  =  T 


1  + 


(7-1) 

2 


(2-44) 


In  high  speed  flows,  a  large  amount  of  kinetic  energy  is  carried  simply  by  the  motion 
of  the  flow.  When  the  flow  decelerates  as  it  approaches  the  wall,  the  kinetic  energy  of  the 
flow  is  converted  to  a  temperature  rise  in  the  fluid. 

The  stagnation  temperature  is  the  temperature  where  all  the  kinetic  energy  is 
converted  to  thermal  energy  in  the  fluid  as  the  fluid  is  slowed  isentropically  from  its  free- 
stream  velocity  to  zero.  However,  in  reality,  there  are  viscous  dissipative  effects  in  the 
boundary  layer  and  sub-layer;  thus,  a  small  amount  of  mass  transport  occurs  normal  to  the 
direction  of  flow.  These  dissipative  effects  cause  the  actual  temperature  to  be  lower  than 
what  is  predicted  by  the  stagnation  temperature;  this  temperature  is  referred  to  as  recovery 
temperature. 

We  can  introduce  a  recovery  factor  r  to  get  the  recovery  temperature,  which  takes 
into  account  the  small  amount  of  mass  transport,  thereby  reducing  the  temperature  slightly 
below  the  stagnation  temperature.  The  recovery  temperature  is: 

Tr  =  T(l  +  r^Y^M2^J  (2-45) 

where  we  would  expect  the  recovery  factor  to  be  less  than  1  in  order  to  slightly  reduce  the 
temperature  below  the  stagnation  temperature.  The  Prandtl  number  is  equal  to  the  viscous 
diffusion  rate  divided  by  the  thermal  diffusion  rate  (White  79): 

ucv  viscous  diffusion  rate  . 

Pr  =  t—P= -  (2-46) 

k  thermal  diffusion  rate 

The  use  of  the  Prandtl  number  is  intuitive  because  the  recovery  temperature  accounts  for 
the  small  amount  of  mass  and  thermal  diffusion  occurring  normal  to  the  flow  at  the 
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stagnation  point.  Schlichting  shows  the  experimentally  found  recovery  factor  for  laminar 
flow  in  air  (Schlichting  335): 

rl  ~  r laminar  ~  V/V  (2-47) 

Schlichting  also  cites  the  recovery  factor  for  turbulent  flow  in  air,  with  experimental  proof 
(Schlichting  713-714): 

rt  ~  r turbulent  ~  V/V  (2-48) 

In  conclusion,  we  can  witness  the  conversion  of  kinetic  energy  in  the  flow  to  a 
temperature  rise  through  some  CFD  results.  Although  the  flow  is  not  isentropic,  this  process 
is  somewhat  reversible.  In  Figures  2-3  and  2-4,  the  incoming  flow  is  Mach  5  and  the  angle  of 
attack  is  45°.  The  flow  decelerates  to  stagnation,  illustrated  by  the  purple  region  near  the 
front  below  the  airfoil's  surface  in  Figure  2-4,  and  creates  a  temperature  rise,  shown  by  the 
red  region  in  Figure  2-5.  The  flow  then  moves  around  the  left  most  tip  of  the  airfoil  and 
accelerates  as  it  reaches  the  top  of  the  airfoil,  shown  by  the  brown  region  above  the  airfoil's 
surface.  This  acceleration  consumes  thermal  energy,  and  thus,  the  flow  cools;  this  is  the 
reversible  process.  This  cool  area  is  the  light  blue  to  dark  blue  region  near  the  center  of  the 
airfoil  above  the  airfoil's  surface  in  Figure  2-5. 
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Figure  2-4  -  Mach  number  contours  for  an  Figure  2-5  -  Temperature  contours  for  an 

airfoil  flying  at  Mach  5  and  45°  angle  of  airfoil  flying  at  Mach  5  and  45°  angle  of 


attack. 


attack. 


Near  the  back  of  the  airfoil,  flow  separation  occurs,  and  a  region  of  recirculation  is  present. 
This  causes  the  flow  to  slow  to  subsonic  speeds,  as  depicted  by  the  purple  region  near  the 
rear  of  the  airfoil  in  Figure  2-4.  This  then  causes  the  fluid  to  again  gain  temperature,  and 
this  can  be  seen  as  the  orange  and  light  red  region  near  the  rear  of  the  airfoil  above  the 
airfoil's  surface. 

Later,  in  Chapter  4,  we  will  see  that  accounting  for  the  kinetic  energy  of  the  flow  by 
using  the  recovery  temperature  instead  of  free  stream  temperature  in  the  formulation  the 
heat  transfer  coefficient  will  reduce  the  heat  transfer  coefficient's  dependence  on  wall 
temperature.  Essentially,  the  recovery  temperature  based  heat  transfer  coefficient  is  weakly 
dependent  on  wall  temperature  when  compared  to  a  standard  formulation  of  heat  transfer 
coefficient,  which  uses  the  free  stream  temperature. 


22 


CHAPTER  3 


LOWER  ORDER  METHODS  TO  CALCULATE  WALL  TEMPERATURE 


Recovery  Temperature 

The  derivation  in  Chapter  2  shows  how  recovery  temperature  arises  from  a 
conversion  of  the  kinetic  energy  in  a  fluid  to  the  temperature  in  a  fluid  when  the  fluid 
decelerates  to  stagnation  from  a  high  speed.  As  witnessed  in  the  course  of  this  research 
project,  the  recovery  temperature,  using  CFD  generated  Prandtl  numbers,  agrees  very  well 
with  steady  state  CFD  results  on  leading  edges.  Therefore,  if  flight  conditions  allow  for 
steady  state,  ideal  gas  assumptions,  the  recovery  temperature  can  be  used  as  a  means  to 
calculate  the  peak  temperature  on  the  leading  edges  of  the  vehicle. 

Steady  state  conditions  rarely  exist  in  real  flights,  however.  Therefore,  the  recovery 
temperature  method  of  calculating  peak  wall  temperature  for  sizing  the  thermal  protection 
system  can  over-predict  the  real  wall  temperature  encountered  in  transient  conditions.  In 
addition,  the  recovery  temperature  only  holds  true  for  ideal  gas  flows,  which  is  not  a  valid 
assumption  in  many  re-entry  problems,  such  as  the  Space  Shuttle  re-entering  the 
atmosphere  at  Mach  20.  The  recovery  temperature  also  fails  to  predict  the  heat  rise  in  the 
fluid  that  is  not  caused  by  kinetic  energy  conversion. 

Over  a  large  portion  of  the  vehicle,  such  as  on  the  upper  and  lower  surfaces  of  the 
wings  and  along  the  fuselage,  the  flow  can  be  thought  of  as  flow  over  a  flat  plate.  In  this  case, 
higher  temperatures  are  a  result  of  viscous  effects  in  the  boundary  layer,  as  opposed  to  the 
stagnation  zone  phenomena  taking  place  on  the  nose  and  wing  leading  edges;  the 
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temperature  rise  in  the  fluid  comes  from  shearing  the  fluid.  The  recovery  temperature  does 
not  account  for  the  shear,  only  the  kinetic  energy,  which  is  signified  by  the  M2  term  in 
Equation  (3-1). 

Tr  =  r(l  +  r^lw2)  (3-1) 

Thus,  using  the  recovery  temperature  cannot  predict  the  amount  of  thermal  protection 
needed  when  the  flow  is  parallel  to  the  surface  of  the  vehicle. 

A  better  method  of  temperature  calculation  over  the  entire  surface  of  the  vehicle 
exists,  which  uses  a  surface  grid  on  the  geometry  and  streamlines  through  the  surface  grid 
to  recognize  when  the  flow  is  perpendicular  or  parallel  to  the  surface  of  the  vehicle. 
Aeroheating  and  Thermal  Analysis  Code  (ATAC)  is  one  of  these  methods. 

Aeroheating  and  Thermal  Analysis  Code 

ATAC  is  widely  used  in  the  aerospace  industry  as  a  method  to  rapidly  predict  wall 
temperatures  on  a  hypersonic  vehicle  geometry.  ATAC  is  higher  order  than  a  simple 
recovery  temperature  calculation,  but  lower  order  than  a  Navier-Stokes  CFD  code. 

ATAC  does  not  grid  the  volume  of  the  solution  domain,  as  a  Navier-Stokes  CFD 
solver  does.  It  uses  a  surface  grid  to  generate  streamlines,  upon  which  it  then  performs 
various  calculations.  Figure  3-1  is  an  example  of  the  ATAC  surface  grid  on  a  reusable 
booster  system  concept  geometry.  As  stated  in  the  ATAC  manual,  the  solution  procedure 
consists  of  an  "inviscid  solution  to  provide  the  boundary  layer  edge  conditions  and  the 
surface  energy  balance  computations...  provide  the  wall  boundary  conditions  for  the 
solution  of  the  integral  equations"  (Strobel,  Tillman  and  King  19).  Because  the  inviscid 
portion  of  the  flow-field  dominates  the  pressure  distribution,  it  can  be  used  to  reflect  the 
state  of  the  boundary  layer  edge  properties.  Subsequently,  the  heat  and  mass  transfer  can 
be  calculated. 
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Figure  3-1  -  ATAC  symmetric  surface  grid  of  the  reusable  booster  concept. 

ATAC  performs  inviscid  streamline  tracing  using  the  following  equation,  known  as 
"Newtonian  approximation  or  the  method  of  steepest  descent"  (Strobel,  Tillman  and  King 

9), 


dr  dr  dw  dr  du 
ds  dw  ds  du  ds 


(3-2) 


Where  t  is  the  streamline  trace,  r  is  the  position  vector  of  some  point  on  the  body.  These 
streamlines  serve  as  the  basis  of  all  future  calculations. 

The  calculation  of  pressure  distribution  along  the  body  can  be  solved  by  any  of  the 
three  choices:  Dahm-Love  pressure  correlations,  Newtonian  pressure,  or  modified 
Newtonian  pressure.  The  Dahm-Love  correlations  are  "an  empirical  extension  and  synthesis 
of  the  modified  Newtonian  correlation,  valid  for  spheres,  but  including  a  correlation  for  flat 
faced  cylinders"  (Strobel,  Tillman  and  King  11).  The  Newtonian  method  for  approximating 
surface  pressure  is  calculated  via  the  following  equations: 


Cp  =  Cp  COS2  T] 


(3-3) 


where  r/  is  the  angle  between  the  wind  vector  and  the  outward  surface  normal.  Cp  is 
calculated  by 


25 


(3-4) 


C*  = 


th 

(IK 


Moo  <  1 


1  + 


- — -M2 

n  1V1  or. 


f  *  = 


© 


Moo  >  1 


(3-5) 


Mi 


Likewise,  the  modified  Newtonian  pressure  is  augmented  with  a  function  of  the  angle 
between  the  wind  vector  and  the  outward  surface  normal,  rj  and  the  Mach  number,  M. 


Cp  =  Cp  cos2  rj  +  F(jj,  M) 


(3-6) 


where 


F(y,  M)  =  [-0.78M-2-27  cos  p  -  o^Se-2'235^00-1)]  sin  r]  (3-7) 


Equation  (3-7)  was  “modified  using  the  empirical  function  of  Andrews  for  low  supersonic 
and  high  subsonic  flight"  (Strobel,  Tillman  and  King  15).  This  relationship  applies  to  sharp 
cone  geometries. 

The  surface  temperature  and  heat  flux  are  calculated  using  the  Momentum/Energy 
Integral  Technique  (MEIT).  The  ATAC  manual  says  the  MEIT  is  “applicable  to  continuum 
flows  and  to  non-separated  boundary  layer  conditions"  (Strobel,  Tillman  and  King  18).  The 
following  three  equations  are  solved  simultaneously  in  the  MEIT  procedure.  The  integral 
momentum  equation, 


rPeUe{Ke  ~  hw)  ds 


(r Peue(ht  e  -  hw)<S> ) 


—  Ch 


JfXy  hw 
ht,e  —  hw 


+ 


(P^)w(ftqw  ~ftw) 
Peue(ht,e  ~ 


(3-8) 


The  integral  energy  equation, 


1  d 
rPeue  ds 

And  the  entrainment  relation, 


(rpeu20)  =  %  + 


(pv)wUjiW  |  HQ  dp 
peu2  peu2  ds 


(3-9) 
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where 

0 

<±  = 
2 

// 

8* 

ht  = 
hw  = 
hr  = 

(D 

Ch  = 

y 

F 

5* 

8 

F 

Fe0  = 


00 


PooUooy2 


=  2  rF/igReg  —  2 


(3-10) 


Boundary  layer  momentum  thickness, 


pu 

Peue 


Friction  factor, 

PeUe 


Wall  shear 


Shape  factor, 

Boundary  layer  displacement  thickness,  J0°°  (l  —  dy 


Total  enthalpy 
Wall  enthalpy 
Recovery  enthalpy 

Boundary  layer  energy  thickness,  J°°  \—  —  )  dy 


Stanton  number, 


_ Q_w _ 

Peue0^r~  ^w) 


Wall  heat  flux 

Shock  radial  location  through  which  the  boundary  layer  edge  streamline  has 
passed 

Entrainment  shape  factor, 

Displacement  thickness  based  on  edge  properties 

Boundary  layer  thickness 

Viscosity 

Momentum  thickness  Reynolds  number,  PeUeQ 


And  the  subscripts  are: 
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w  =  Wall 

e  =  Boundary  layer  edge 

i,w  =  Properties  obtained  through  an  isentropic  expansion  from  the  stagnation  point 
to  the  location  pressure 

The  following  excerpt  from  the  ATAC  manual  describes  each  equation's  purpose  towards 
the  solution: 

The  momentum  equation  solution  dictates  the  boundary  layer  parameters  upon 
which  the  following  phenomena  are  based:  transition  onset  and  location, 
transitional  intermittency,  surface  roughness  effects,  and  turbulent  boundary  layer 
shape  factors.  The  energy  equation  solution  dictates  the  convective  heat  transfer 
subject  to  the  four  parameters  above.  The  boundary  layer  parameter,  which  dictates 
convective  mass  transfer,  is  not  evaluated.  Rather,  mass  transfer  is  determined  by 
analogy  to  the  energy  boundary  layer  solution.  The  entrainment  relation,  Equation 
(3-10),  provides  a  means  of  determining  the  boundary  layer  properties,  which  are 
essential  boundary  conditions  for  the  solution  of  Equations  (3-8)  and  (3-9).  The 
boundary  layer  edge  thermodynamic  state  is  determined  by  lookup  on  pressure  and 
entropy  in  a  real-gas  Mollier  table.  Pressure  is  known  from  the  inviscid  flow 
solution,  and  entropy  is  calculated  from  consideration  of  the  bow  shock  shape  and 
boundary  layer  mass  entrainment.  (Strobel,  Tillman  and  King  19) 

These  empirical  correlations  that  compute  pressure  break  down  at  high  angles  of 

attack  and  large  amounts  of  separation,  such  as  on  sharp  leeward  angles.  Due  to  the  rocket- 

back  maneuver  and  the  high  angles  of  attack  during  re-entry,  these  empirical  relationships 

become  questionable  during  a  large  portion  of  the  flight  regime.  According  to  the  ATAC 

manual,  “since  the  methodology  does  not  allow  for  separation  or  flow  across  streamlines, 

the  accuracy  at  these  large  angles  [of  attack]  are  questionable,  especially  on  the  leeward 

side"  (Strobel,  Tillman  and  King  9). 

A  method  capable  of  calculating  at  high  angles  of  attack  is  to  use  a  Navier-Stokes 
based  CFD  code  that  has  a  means  of  directly  modeling  eddies  and  a  method  of  calculating 
turbulence  transition  location  and  intensity.  As  stated  previously,  the  ultimate  method  is 
transient,  fully  coupled  CFD,  but  this  method  requires  more  time  than  is  feasible  for  any 
practical  problem.  Thus,  we  employ  steady  state  Navier-Stokes  CFD. 
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Steady-State  Navier-Stokes  Computational  Fluid  Dynamics 

If  cost  and  time  were  not  relevant  constraints  in  an  engineering  problem  then  by  far, 
the  best  way  to  solve  the  transient  heat  transfer  problem  for  a  thermal  protection  system 
would  be  using  a  fully  coupled,  transient,  conjugate  heat  transfer  analysis.  Several  problems 
arise  with  this  utopian  approach.  First,  the  grid  size  is  extremely  large  because  the  internal 
geometry  must  now  be  modeled,  including  the  thickness  of  all  solids  as  well  as  the  fluid  in 
the  internal  bays.  Secondly,  the  time-steps  must  be  extremely  small,  on  the  order  of  10-6 
seconds,  for  the  stability  of  the  transient  Navier-Stokes  Equations.  Third,  the  wall  clock  time 
required  for  such  small  time  steps  over  a  seven  minute  flight  is  immense.  The  Department 
of  Defense's  (DoD)  high  performance  computing  center  has  a  maximum  allowable  wall 
clock  time  of  two  weeks,  which  is  insufficient  for  the  TPS  problem  at  hand.  To  get  around 
the  DoD  wall  clock  time  limitation,  an  attempt  was  made  to  run  a  sequence  of  transient 
analyses  by  restarting  the  computation  on  a  test  model  every  two  weeks.  Restarting  the 
transient  problem  proved  to  cause  numerical  divergence,  rendering  this  “fix"  useless. 
Coupling  the  above  problems,  we  can  see  that  performing  a  transient,  conjugate  heat 
transfer  analysis  is  not  possible  with  current  computer  technology.  Thus,  we  seek  a  set  of 
steady-state  solutions  generated  by  a  Navier-Stokes  CFD  code  as  an  effective  engineering 
solution  to  the  real  world  computational  limitations. 

Steady  state  Navier-Stokes  CFD,  by  definition,  does  not  include  transient  effects  in 
the  flow-field.  These  transient  phenomena  may  include  vortex  shedding,  an  undeveloped 
boundary  layer,  or  movement  of  the  shock  layer,  and  can  have  an  effect  on  the  temperature 
distribution  along  the  wall  of  the  vehicle.  The  surface  temperatures  of  the  vehicle  are 
dominated  by  both  the  conversion  of  the  flow's  kinetic  energy  to  a  temperature  rise  at  the 
leading  edges  and  the  viscous  shearing  of  the  fluid  along  flow-parallel  surfaces,  not 
necessarily  by  transient  effects. 
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It  is  quite  common  in  the  CFD  community  to  use  steady  state  or  quasi-steady  state 
analyses,  it  is  based  on  the  fact  that  “the  flow-field  adapts  at  a  very  fast  rate  to  changes  in 
the  free  stream  conditions  as  the  vehicle  is  accelerating.  This  is  in  contrast  to  the  solid  side 
temperature  profile  which  gradually  changes  throughout  the  flight,  never  reaching  steady 
state"  (Fricker,  Mendoza  and  Catton  4). 

Commercial  CFD,  with  readily  available  computational  resources,  is  capable  of 
carrying  out  the  steady  state  analysis.  Ansys  Fluent  uses  the  full  set  of  Navier-Stokes 
Equations,  which  are  comprised  of  the  Continuity,  Momentum,  and  Energy  Equations.  The 
continuity  equation  is  solved  in  the  following  form  (Ansys,  Inc.  27): 

^  +  V-(pv)=5m  (3-11) 

The  source  mass  Sm  is  0  for  this  application.  The  momentum  equation  is  given  in  general 
form  (Ansys,  Inc.  28): 

—  (pv)  +  V  ■  (pvv)  =  -Vp  +  V  ■  f  +  pg  +  F  (3-12) 

ot 

where  p  is  the  static  pressure,  f  is  the  stress  tensor,  pg  is  the  gravitational  body  force,  and  F 
is  the  external  body  force.  The  last  two  terms  in  Equation  (3-12)  are  negligible.  “The  stress 
tensor  is  given  by: 

r  2  l 

f  —  ii  (Vv  +  Vvr)  —  -  V  ■  v/  (3-13) 

where  p  is  the  viscosity,  /  is  the  unit  tensor,  and  Vv  +  Vvr  is  the  effect  of  volume  dilation" 
(Ansys,  Inc.  28).  The  energy  equation  is  modeled  in  Ansys  Fluent  “using  the  concept  of  the 
Reynolds'  analogy  to  turbulent  momentum  transfer.  The  'modeled'  energy  equation  is  as 
follows"  (Ansys,  Inc.  Ill): 

d  d  Nn  d  \f  CV,pt\  dT  s  \  ] 

Tt(pE)  +Wlut(pE  +  P)]  =  Wj  [(*  +~^)^  +  “I (*</),„ J  +  (3-U) 
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where  Sh  is  the  source  term  and  the  subscript  “eff"  means  effective.  The  thermal 
conductivity  is  k,  E  is  the  total  energy,  and  (Tiy)eff  is  the  deviatoric  stress  tensor,  given  by 
(Ansys,  Inc.  Ill): 

_  (duj  duA  2  duk 

OoOeff  -  ^eff  (^x.  +  dxJ  3  Meff  g^  $ij  (3'15) 

Fluent  has  a  variety  of  turbulence  models,  from  the  one-equation  Spalart-Allmaras 
to  the  seven-equation  Reynold's  Stress  model.  For  this  application,  the  Transition  SST 
model  was  chosen.  Only  two  transition  models  are  available  in  Ansys  Fluent:  Transition  SST 
and  Reynolds  Stress.  The  Reynolds  Stress  model  has  seven  equations,  more  than  what  is 
being  solve  for  the  general  flow-field.  In  order  to  produce  more  rapid  results,  the  transition 
model  with  fewer  equations  was  chosen. 

This  Transition  SST  model  is  “based  on  the  coupling  of  the  SST  k  —  co  transport 
equations  with  two  other  transport  equations,  one  for  the  intermittency  and  one  for  the 
transition  onset  criteria,  in  terms  of  momentum-thickness  Reynolds  number"  (Ansys,  Inc. 
142).  According  to  White,  “there  is  no  fundamental  theory  of  transition,  but  there  are 
experiments  and  correlations  which  try  to  predict  the  final  onset  of  fully  turbulent  flow, 
such  as  Rex  tr  or  Reetrt  as  a  function  of  the  following  parameters:  pressure  gradient,  free- 
stream  turbulence,  wall  roughness,  mach  number,  wall  suction  or  blowing,  and  wall  heating 
or  cooling"  (White  378).  The  model  used  by  Ansys  uses  an  empirical  correlation  by  Langtry 
and  Menter  “which  covers  standard  bypass  transition  as  well  as  flows  in  low  free-stream 
turbulence  environments"  (Ansys,  Inc.  142). 

The  solution  method  employed  in  this  thesis  uses  steady  state  CFD  solutions  to 
predict  heat  transfer  coefficients  that  are  shown  to  be  only  weakly  dependent  on  wall 
temperature.  Thus,  transient  effects  are  neglected.  Standard  calculation  of  heat  transfer 
coefficients,  in  the  form  of  the  following  equation, 
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(3-16) 


require  knowledge  of  the  actual  wall  temperature  Tw  in  transient  flight,  as  well  as  data 
extraction  from  the  edge  of  the  boundary  layer  for  the  temperature  at  that  location. 
Unfortunately,  the  wall  temperature  is  unknown  a  priori,  as  we  are  not  performing  a 
transient  calculation.  In  addition,  in  regions  of  highly  separated,  re-circulating  flow,  such  as 
the  back  of  the  vehicle  when  the  vehicle  flies  at  angles  of  attack  up  to  90°,  there  is  a  question 
as  to  where  the  edge  of  the  boundary  layer  data  should  be  extracted.  An  alternate  method  to 
predicting  standard  heat  transfer  coefficients  is  to  use  one  based  on  recovery  temperature 

The  recovery  temperature  is  shown  in  Chapter  4  to  be  only  weakly  dependent  on 
wall  temperature.  Therefore,  a  simple  steady  state  analysis  can  be  generated,  the  recovery 
temperature  calculated,  and  a  heat  transfer  coefficient  then  generated  for  use  in  subsequent 
analyses.  The  workflow  follows  this  sequence: 

1.  Generate  steady  state  CFD  results  at  each  time-step  using  a  re-radiating, 
adiabatic  wall,  also  known  as  radiation  equilibrium  wall  temperature  (REWT). 

2.  Export  the  results  in  Tecplot  format.  This  format  is  text  based  lists  the  x,  y,  and  z 
coordinates,  then  any  desired  data  in  the  additional  columns.  The  data  being 
exported  is: 

a.  Far-field  temperature 

b.  Far-field  mach  number 

c.  Wall  temperature 

d.  Heat  flux 

e.  Effective  Prandtl  number.  This  number  includes  the  turbulence  effects 
and  Fluent  calculates  this  number  by  the  effective  thermal  conductivity 
and  effective  viscosity,  which  also  account  for  the  turbulent  effects. 
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3.  Use  a  Matlab  code  specifically  written  for  this  procedure  to  import  the  Tecplot 
results  and  calculate  the  recovery  temperature  at  each  point.  This  code  uses  the 
heat  flux,  wall  temperature,  and  Prandtl  number  at  each  element  on  the  wall. 

4.  Export  the  results  in  a  comma  delimited  format  for  use  with  the  ABAQUS  FEA 
heat  transfer  analysis 

Before  the  applied  heat  transfer  analysis  can  be  performed  on  the  full  three- 
dimensional  geometry,  the  dependency  of  the  heat  transfer  coefficient  with  respect  to  wall 
temperature  needs  to  be  quantified.  This  dependency  needs  to  be  quantified  in  order  to 
validate  the  decoupled  approach  to  the  solution  of  the  transient  TPS  temperature  profile. 
For  this  process  of  quantification,  a  simplified,  two-dimensional  case  is  chosen  for  reduced 
computation  time  and  less  complexity  in  code. 
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CHAPTER  4 


QUANTIFYING  THE  DEPENDENCY  OF  THE  HEAT  TRANSFER  COEFFICIENT  WITH  RESPECT 

TO  WALL  TEMPERATURE 


Weakly  Dependent  Heat  Transfer  Coefficient 

The  recovery  temperature  based  heat  transfer  coefficient  is  only  weakly  dependent 
on  wall  temperature  because  it  accounts  for  the  kinetic  energy  of  the  incoming  flow.  Thus, 
having  calculated  the  heat  transfer  coefficient  for  one  case,  the  wall  temperature  need  not 
be  known  for  the  other  cases;  furthermore,  the  heat  transfer  coefficient  can  then  be  applied 
to  a  transient  heat  transfer  analysis  on  a  thermal  protection  system  and  the  transient 
temperature  profile  can  be  produced  with  a  higher  degree  of  confidence.  The  data 
presented  in  this  chapter  will  show  that  the  recovery  temperature  based  heat  transfer 
coefficient  is  weakly  dependent  on  wall  temperature  for  a  wide  range  of  flight  conditions. 

The  Wing  Model 

A  simple  two-dimensional  model  of  a  wing  provides  the  geometry  for  the  variability 
quantification  analysis.  Figure  1-1  depicts  the  geometry,  which  is  a  NACA  23012  airfoil.  The 
two-dimensional  model  was  chosen  in  order  to  simplify  any  quantification  of  the  variation 
in  heat  transfer  coefficient  with  respect  to  wall  temperature  and  to  neglect  any  three- 
dimensional  flow  dissimilarities.  Additionally,  solution  time  is  massively  reduced,  as  the 
grid  size  is  cut  by  several  orders  of  magnitude:  the  three-dimensional  grid  used  later  in  this 
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research  has  26,245,437  elements  and  6,083,498  nodes  while  the  two-dimensional  grid  of 
the  wing  has  only  44,550  elements  and  44,880  nodes. 

The  grid  spacing  in  the  direction  normal  to  the  surface  had  to  be  created  small 
enough  near  the  wall  in  order  to  capture  the  velocity  gradient  in  the  boundary  layer.  The 
grid  spacing  is  then  expanded  exponentially  to  the  far-field  domain,  which  is  20  chord- 
lengths  away.  The  expansion  in  grid  spacing  serves  as  a  way  to  reduce  mesh  size  in  areas 
that  steep  gradients  do  not  exist,  such  as  beyond  the  boundary  layer.  Figure  4-1  shows  the 
mesh  near  the  wall. 


Figure  4-1  -  Grid  of  the  NACA  23012  airfoil  used  in  the  two-dimensional  wall  temperature 
dependency  analysis. 

The  grid  statistics  are  as  follows.  The  maximum  aspect  ratio  is  161.6,  which  is  high, 
but  the  short  side  of  the  cell  is  parallel  to  the  direction  of  high  temperature  and  velocity 
gradients  in  the  boundary  layer  and  the  gradients  along  the  long  side  are  much  less  severe. 
Figure  4-2  shows  an  enlarged  view  of  the  area  of  high  aspect  ratio.  The  high  aspect  ratio  is  a 
result  of  the  expansion  layer  used  to  resolve  the  velocity  and  temperature  profiles  in  the 
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boundary  layer.  In  order  to  reduce  the  aspect  ratio,  a  large  number  of  divisions  in  the 
direction  parallel  to  the  surface  of  the  body  would  be  necessary.  This  would  greatly  increase 
the  cell  count  and  therefore  increase  the  time-to-solution. 


Figure  4-2  -  Enlarged  view  of  high  aspect  ratio  region. 


The  maximum  cell  skewness  is  85.3%,  well  within  the  limit  for  a  double  precision  solver. 
Ansys  calculates  skewness  via  the  following  equation: 

n1  f@max  ~  @e  @e  ~  @min\ 

Skewness  =  max  — — — — , - - - 

V 180  -9e  9e  ) 

where  9e  is  the  equiangular  face  or  cell,  which  is  60°  for  tetrahedrons  and  triangles  and  90° 
for  quadrilaterals  and  hexahedrons.  Figure  4-3  depicts  9max  and  9min. 

@max 

cell 

Figure  4-3  -  Graphic  of  9max  and  9min. 
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The  maximum  skewness  occurs  near  the  trailing  edge  of  the  airfoil  and  exists  in  the 
boundary  layer.  Because  the  grid  is  structured,  the  skewness  is  a  result  of  the  steep  change 
in  curvature  of  the  trailing  edge  of  the  airfoil  versus  the  gradual  change  in  curvature  of  the 
opposing  boundary. 


Figure  4-4  -  Enlarged  view  of  the  region  with  the  highest  skewness. 

The  full  set  of  Navier-Stokes  equations  is  modeled;  these  equations  appear  in  the 
previous  chapter  as  Equations  (3-11)  through  (3-15).  The  turbulence  model,  Transition 
SST,  is  discussed  in  Chapter  3.  It  is  a  four-equation  turbulence  model  that  uses  the  k  —  co 
transport  equations  and  an  Ansys  proprietary  transition  model  derived  from  empirical 
curve  fits.  Air  is  modeled  as  an  ideal  gas.  The  specific  heat  cp  is  described  by  a  piecewise- 
polynomial  with  two  ranges  and  eight  coefficients,  a  built-in  Fluent  model.  Table  4-1 
summarizes  the  piecewise-polynomial  coefficients  for  each  range.  The  model  for  thermal 
conductivity  and  viscosity  is  kinetic  theory. 
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Table  4-1  -  Ranges  and  coefficient  values  for  the  piecewise  polynomial  model  of  specific 
heat  Cp. 


Range  1:  100  K  <  T  <  1000  K 

Coefficient  # 

1 

2 

3 

4 

5 

6 

7 

8 

Value 

1161 

-2.369 

0.01486 

-5.035E-5 

9.929E-8 

-1.11E-10 

6.54E-14 

-1.57E-17 

Range  2:  1000  K  <T  <  3000  , 

K 

Coefficient  # 

1 

2 

3 

4 

5 

6 

7 

8 

Value 

-7070 

33.71 

-0.058 

5.422E-5 

-2.337E-8 

9.24E-12 

-1.57E-15 

1.11E-19 

There  are  several  boundary  conditions  that  were  chosen  to  bound  the  weak 
dependence  problem.  For  Study  A,  Table  4-2  lists  the  various  boundary  conditions.  These 
boundary  conditions  were  chosen  to  represent  the  flight  trajectory  of  the  applied  problem 
described  in  Chapter  5.  The  boundary  conditions  listed  in  Table  4-1  compare  isothermal 
walls  with  walls  that  have  a  temperature  gradient:  the  REWT  wall,  the  Ti-6-4  wall,  and  the 
Inconel-617  wall.  The  Ti-6-4  wall  is  the  titanium  6-aluminum  4-vanadium  alloy  with  a  wall 

/c  Q 

thickness  of  0.125  in.  Ti-6-4  properties  are  a  density,  p,  of  4430—,  a  specific  heat  capacity, 
Cp,  of  526  kJ  — ,  and  a  thermal  conductivity,  k,  of  6.7^^  (Matweb,  LLC.).  The  Inconel-617 
wall  is  the  Special  Metals®  Alloy  617  with  a  wall  thickness  of  0.125  in.  The  material 

kci  J 

properties  are  a  density,  p,  of  8360^-,  a  specific  heat  capacity,  Cp,  of  419  k  R,  and  a 

thermal  conductivity,  k,  of  13.6^^;  (Matweb,  LLC.).  Both  of  these  metals  are  typical 
structural  alloys  used  in  high  speed  aircraft. 
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Table  4-2  -  Boundary  conditions  for  CFD  jobs  used  in  Study  A  of  quantifying  the  variation  in 
the  weak  dependence  assumption. 


BOUNDARY  CONDITIONS  -  STUDY  A 

Angle  of  Attack,  a  (°) 

Wall 

Emissivity 

Mach  Number 

0 

200  K 

0.8 

2 

15 

600  K 

0.9 

5 

45 

1000  K 

90 

REWT 

Ti-6-4 

Inconel-617 

The  boundary  conditions  listed  in  Table  4-2  are  combined  to  produce  72  unique 
CFD  runs  for  Study  A,  after  removing  duplicate  boundary  conditions.  Table  B-l  in  Appendix 
B  contains  a  table  of  all  the  various  boundary  conditions  for  each  job  code.  In  addition  to 
these  boundary  condition,  a  constant  altitude  of  30  km  is  used,  which  results  in 
Poo  =  1172  Pa  and  T =  226.65  K  using  the  1976  U.S.  Standard  Atmosphere. 

In  addition  to  Study  A,  which  compares  isothermal  walls  to  radiation  adiabatic 
walls,  there  is  a  separate  study,  Study  B,  which  measures  the  dependency  of  the  heat 
transfer  coefficient  with  respect  to  wall  temperature,  using  isothermal  wall  boundary 
conditions  only.  In  this  case,  the  boundary  conditions  are  varied  only  at  the  wall.  There  are 
no  changes  to  Mach  number  or  angle  of  attack.  The  same  altitude  of  30  km  is  used,  again 
resulting  in  =  1172  Pa  and  =  226.65  K  using  the  1976  U.S.  Standard  Atmosphere. 
The  Mach  number  is  5  and  the  thermal  boundary  conditions  at  the  wall  were  varied  from 
200  K  to  1000  K  in  100  K  increments,  resulting  in  9  unique  analyses. 

For  all  the  analyses  reported  in  this  technical  document,  the  pressure  far-field 
turbulent  intensity  is  set  to  0.01%  and  the  turbulent  viscosity  ratio  is  set  to  10.  Roy  and 
Blottner  investigated  several  turbulence  models,  including  a  Wilcox  k  —  go  model,  which  is 
the  basis  for  the  Ansys  k  —  go  model  (Ansys,  Inc.  126).  Part  of  Roy  and  Blottner's 
investigation  included  a  sensitivity  analysis  of  Wilcox's  k  —  go  turbulence  model  to  free- 
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q  (W/m2) 


stream  turbulence  levels  as  well  as  free-stream  turbulence  viscosity  ratios  on  a  Mach  8 
sharp  cone.  Figure  4-5  shows  the  heat  flux  as  a  function  of  distance  along  the  cone  for  each 
case.  Table  4-3  summarizes  the  variations  in  viscosity  ratio  and  level  of  intensity  in  all  the 
different  cases. 


Figure  4-5  -  Sensitivity  of  the  Wilcox  k  —  go  turbulence  model  to  the  free-stream  turbulence 
levels  for  the  Mach  8  sharp  cone  (Roy  and  Blottner  321). 
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Table  4-3  -  Summary  of  the  cases  studied  in  Figure  4-5,  showing  the  variations  in  free- 
stream  turbulent  viscosity  ratio  and  turbulence  intensity.  Y  denotes  a  case  that  was  run 
(Roy  and  Blottner  321). 


Tu>  % 

A 

B 

C 

D 

E 

Case 

0.001 

0.01 

0.1 

1.0 

10.0 

i 

1  x  10-5 

2 

0.001 

y 

Y 

3 

0.1 

Y 

Y 

Y 

4 

1.0 

Y 

Y 

Y 

5 

10.0 

Y 

Y 

Y 

Y 

The  solution  methods  are  comprised  of  an  implicit  formulation  and  an  Advection 
Upstream  Splitting  Method  (AUSM)  flux  type.  The  AUSM  flux  type  is  more  robust  for 
problems  with  shock  waves:  according  to  Blazek,  the  “advantage  can  be  seen  in  only  a 
moderately  increased  numerical  effort  but  a  much  better  resolution  of  shocks  and  boundary 
layers,  as  compared  to  the  central  scheme  with  scalar  artificial  dissipation"  (Blazek  101). 
Blazek  also  states  that  “AUSM  proved  to  deliver  a  crisp  resolution  of  strong  shocks  and 
accurate  results  for  boundary  layers"  (Blazek  105).  The  other  option  for  flux  type  is  Roe's 
Forward  Difference  Scheme,  which  has  disadvantages  with  chemically  reacting  flows 
(Blazek  110).  Although  in  this  study  the  flow  never  gets  hot  enough  to  chemically  react, 
AUSM  was  chosen  for  its  greater  range. 

The  method  of  discretization  and  evaluation  of  the  gradient  is  Green-Gauss  cell 
based  gradient  computation.  This  is  used  because  it  is  slightly  faster  than  a  least-squares 
approach:  for  Green-Gauss  cell  based  gradient  calculation  “no  additional  data  structures  are 
needed  for  the  reconstruction  of  gradients"  (Blazek  164).  The  flow  discretization  is  second 
order  upwind,  but  the  turbulent  kinetic  energy,  specific  dissipation  rate,  and  intermittency 
are  all  discretized  by  the  first  order  upwind  method. 

The  solutions  were  judged  as  converged  when  the  integrated  total  surface  heat  rate 
stopped  changing  after  three  significant  figures  from  one  iteration  to  the  next.  It  was 
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observed  that  the  integrated  total  heat  rate  was  the  last  metric  to  converge;  often,  mass, 
momentum,  and  energy  had  converged  to  six  significant  figures  before  the  heat  flux 
converged  to  three  significant  figures.  Figure  4-6  provides  as  a  sample  graph  of  the 
residuals  versus  iteration  number. 
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Figure  4-6  -  Example  of  residuals  versus  iteration  number.  “Intermit"  and  "retheta"  are 
turbulence  parameters  for  turbulence  intermittency  and  Ree. 


Grid  Convergence  Study 

A  grid  convergence  study  is  necessary  to  ensure  that  the  results  are  independent  of 
the  mesh.  Examples  of  parameters  that  are  highly  dependent  on  the  mesh  are  shock 
location  and  structure,  resolution  of  the  boundary  layer,  and  turbulence  intensity.  A  grid 
convergence  study  is  conducted  on  the  standard  mesh,  a  once  refined  mesh,  and  a  twice 
refined  mesh.  The  once  refined  mesh  cuts  each  quad  of  the  standard  mesh  into  four  quads 
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and  the  twice  refined  mesh  cuts  each  quad  on  the  once  refined  mesh  into  16  quads.  Table  4- 


4  summarizes  the  meshes  and  their  sizes. 


Table  4-4  -  Summary  of  the  grids  used  in  the  grid  convergence  study. 


Mesh 

Nodes 

Elements 

Nodes  Along  the 
Wall 

Standard 

44,880 

44,550 

330 

Once  Refined 

178,860 

178,200 

660 

Twice  Refined 

714,120 

712,800 

1,320 

The  standard  grid  proved  to  have  sufficient  resolution;  Figure  4-7  shows  the  wall 
temperature  as  a  function  of  path  length,  which  is  distance  along  the  surface  of  the  wall,  and 
focuses  on  the  region  near  the  leading  edge,  where  the  peak  temperature  occurs.  The 
leading  edge  is  a  region  of  interest.  This  graphic  shows  that  the  answer  changes  only  a  small 
quantity  when  the  grid  is  refined  by  a  factor  of  four.  For  example,  the  predicted  peak 
temperature  on  the  coarsest  grid  is  approximately  1.5%  lower  than  the  predicted  peak 
temperature  of  the  finest  grid.  Figure  4-8  shows  the  temperature  profile  of  the  entire  wing. 
In  this  plot,  we  observe  that  all  three  temperature  curves  along  the  wall  of  the  vehicle  look 
nearly  identical,  except  at  the  trailing  edge. 
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Figure  4-7  -  Wall  temperature  as  a  function 
of  path  length  for  each  grid,  zoomed  on  the 
leading  edge. 


Figure  4-8  -  Wall  temperature  as  a  function 
of  path  length  for  each  grid  for  the  entire 
surface  of  the  airfoil. 


There  is  some  discrepancy  in  predicted  temperatures  at  the  trailing  edge, 
approximately  2.75 m  path  length.  This  is  due  to  a  small  amount  of  turbulent  eddying 
present,  as  shown  in  Figure  4-9,  as  well  as  the  expansion  wave  that  occurs.  Because  the  heat 
fluxes  on  this  part  of  the  vehicle  are  low,  the  answer  in  this  region  is  not  significant  from  a 
TPS  or  wing  design  perspective. 
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Figure  4-9  -  Enlarged  view  of  the  results  near  the  trailing  edge  of  the  airfoil,  depicting  the 
recirculation  zone.  Note  that  the  dual  vortex  pattern  is  seen  only  in  a  steady  calculation:  a 
transient  calculation  would  show  the  vortices  alternating  and  flowing  downstream. 


Quantifying  the  Wall  Temperature  Dependency  in  the  Heat  Transfer  Coefficient 

In  study  B,  the  nine  analyses  were  created  to  study  the  effect  of  wall  temperature  on 
recovery  temperature  based  heat  transfer  coefficient  alone.  These  nine  analyses  consisted 
of  isothermal  boundary  conditions  at  the  wall,  starting  from  200  Kt  in  100  ^increments,  up 
to  1000  K  The  results  were  exported  in  a  delimited  format  from  Fluent  and  brought  into  a 
simple  Excel  spreadsheet.  Recovery  temperature  was  calculated  by  Equation  (2-45)  and  the 
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recovery  factor  was  calculated  by  Equation  (2-48).  The  Prandtl  number  used  was  the 
effective  Prandtl  number,  which  combines  both  the  laminar  and  turbulent  Prandtl  numbers. 

The  data  had  to  be  plotted  in  a  useful  order.  For  example,  plotting  wall  temperature 
versus  x-location  on  the  grid  is  senseless  because  the  plot  would  reflect  the  upper  surface, 
then  the  lower  surface,  then  the  upper,  or  some  combination  thereof.  Essentially,  plotting 
data  vs.  x-location  does  not  discriminate  what  surface  of  the  airfoil  the  data  lays  (such  as  the 
upper  or  lower  surface).  The  following  method  was  used  to  prepare  the  data  for  plotting: 

1.  Import  the  data  into  an  Excel  spreadsheet 

2.  Find  the  geometric  center  of  the  grid  by  using  the  following  equation: 

_  max(x)  +  min(x) 

xloc  ~  o 


max(y)  +  min(y) 

yioc  = - 2 -  (4"2) 

where  x  and  y  are  the  x  and  y  coordinates  of  the  nodes  along  the  wall. 

3.  Calculate  the  angle  of  each  node  with  the  positive  horizontal  by  the  following 
equation: 

6i  =  atan2(xj  -  xloc,yi  -  yioc)  (4-3) 

where  i  goes  from  1  to  the  number  of  nodes  along  the  wall  and  atan2is  the 
arctangent  of  the  specific  coordinate  between  -  n  and  n. 

4.  Sort  the  data  by  the  angle,  from  smallest  to  largest 

5.  Calculate  the  distance  along  path,  s,  by  this  equation: 

Si  =  V (x;  -  Xj.i)2  +  (y*  -  y£_i)2  (4-4) 

at  i  =  1,  Si  =  0. 

6.  Calculate  the  total  path  /?* : 


Pi  ~  si  +  Pi-1 


(4-5) 
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at  i  =  1,  p1  =  0.  The  total  path  signifies  the  path  length  along  the  surface  of  the 
airfoil. 

7.  Plot  the  data  versus  total  path. 

For  Study  B,  a  metric  was  chosen  to  demonstrate  the  variability  of  the  heat  transfer 
coefficient  as  a  function  of  wall  temperature  and  nondimensionalized  heat  flux.  At  each 
element  along  the  wall,  the  heat  flux  is  nondimensionalized  by  calculating  it's  heat  flux 

versus  the  element  with  the  highest  heat  flux  along  the  wall.  This  is  simply  — —  for  each 

Qmax 

element.  This  result  is  plotted  against  variation  in  the  heat  transfer  coefficient,  ehr  from  the 
mean.  For  this  case,  the  mean  is  calculated  as: 

9 

K  =  |£  hn  (4'6) 

i=  1 

The  variation  is  also  calculated  for  each  element  along  the  wall.  For  example,  at  M  =  5  and 
a  —  0,  let  us  say  that  the  forward  most  element,  the  one  at  the  leading  edge,  has  the  highest 

heat  flux  of  all  the  elements  along  the  wall.  This  element's  — —  =  1  and  it  is  compared  to 

Qmax 

each  leading  edge  element  for  each  wall  boundary  condition,  from  200  K  <TW  <  1000  K. 
Figure  4-10  illustrates  the  results  of  Study  B. 

The  results  of  Study  B  suggest  a  promising  outlook  for  the  use  of  a  recovery 
temperature  based  heat  transfer  coefficient.  We  see  here  that  in  all  cases  except  for  the 
1000  K  case,  the  elements  above  5%  of  the  maximum  heat  flux,  or  0.05  <  — —  <  1,  the 

Qmax 

variation  in  heat  transfer  coefficient  is  between  -10%  and  10%,  or  -10%  <  shr  <  10%.  At 
the  highest  heat  flux  element,  the  band  in  variation  is  much  tighter:  —3%  <  £hr  <  2%  for 

— —  =  1.  This  means  that  for  areas  of  concern  to  the  TPS  engineer,  where  the  heat  flux  is 

Qmax 

high,  the  heat  transfer  coefficient  is  fairly  independent  of  wall  temperature.  This  result  is 
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significant;  without  it,  the  computational  method  used  in  this  research  would  be  of  minimal 


value  to  the  TPS  design  engineer. 


Variation  in  hr  (%) 


Figure  4-10  -  Comparison 


of  — — —  vs.  variation  in  hr  for  each  isothermal  boundary  condition. 

Qmax 


Clearly,  the  heat  transfer  coefficient  become  less  trustworthy  at  the  trailing  edge  of 
the  vehicle,  where  a  band  in  variation  is  observed  to  be  —60%  <  £hr  <  120%  but  this  is  for 

0.003  <  — —  <  0.032,  elements  that  contribute  almost  no  heat  flux  to  the  wall,  which  is 

Qmax 

illustrated  in  Figure  4-10.  This  larger  difference  is  a  result  of  the  flow  above  and  below  the 
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airfoil  coming  back  together,  turbulent  effects  in  the  wake,  and  expansion  waves;  the  flow 
physics  are  much  more  complex  in  this  region. 

A  much  different  result  in  wall  temperature  dependence  of  the  heat  transfer 
coefficient  is  observed  when  a  thermal  gradient  in  the  wall  is  introduced.  Figure  4-10  is  re¬ 
plotted  to  include  a  Radiation  Equilibrium  Wall  Temperature  (REWT)  wall.  It  is  profoundly 
evident  that  the  heat  transfer  coefficient  becomes  sensitive  to  the  thermal  gradient  in  the 
wall,  nearly  linearly  sensitive  for  — —  >  0.3,  but  still  a  good  prediction  for  — —  =  1,  where 

Qmax  Qmax 

the  variation  in  heat  transfer  coefficient  is  ±  2%,  as  seen  in  Figure  4-11. 


Variation  in  hr  (%) 


Figure  4-11  -  Comparison  of  — —  vs.  variation  in  hr  for  each  isothermal  boundary  condition 

Qmax 

and  one  REWT  boundary  condition. 
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The  72  steady  state  analyses  in  Study  A  produced  similar  results  to  Study  B.  The 
walls  with  finite  thickness  and  a  material  property,  both  Inconel-617  and  Ti-6-4,  appeared 
to  have  similar  thermal  characteristics  as  the  REWT  wall.  That  is,  they  had  similar  peak 
temperatures  and  thermal  gradients.  The  variation  in  heat  transfer  coefficient  as  a  function 
of  path  length  along  the  surface  of  the  wing  have  comparable  curves  for  the  similar 
boundary  condition  types  for  the  majority  of  light  conditions,  except  for  the  90°  case.  The 
following  set  of  figures,  Figure  4-12  through  Figure  4-16,  show  the  percentage  variation  in 
heat  transfer  coefficient  as  a  function  of  path  length  along  the  wall  for  some  selected  cases. 
They  are  grouped  by  Mach  number  and  angle  of  attack;  the  remaining  sets  of  plots  which 
investigate  the  percentage  variation  in  heat  transfer  coefficient  versus  path  length  for  Study 
A  are  illustrated  in  Appendix  A. 


Figure  4-12  -Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations 
in  wall  boundary  conditions  at  Mach  2  and 
a  =  0°. 


Figure  4-13  -  Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations 
in  wall  boundary  conditions  at  Mach  5  and 
a  =  0°. 
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Alpha  =  15,  Mach  =  5  Alpha  =  45,  Mach  =  5 


Figure  4-14  -  Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations 
in  wall  boundary  conditions  at  Mach  5  and 
a  =  15°. 


Figure  4-15-  Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations 
in  wall  boundary  conditions  at  Mach  5  and 
a  =  45°. 


As  seen  in  Figures  4-12  through  4-15,  we  see  that  the  REWT  boundary  condition 
closely  approximates  the  heat  transfer  coefficient  that  the  conducting  walls  have.  We  can 
also  see  that,  for  the  Mach  2  cases,  the  heat  transfer  coefficient  has  a  large  amount  of 
dependence  on  wall  temperature  at  the  leading  edge,  a  path  length  of  0  m,  which  is  seen  in 
Figure  4-12.  This  dependence  is  due  to  the  calculation  of  the  heat  transfer  coefficient: 


hr 


where  Tw  ~  Tr  for  the  Mach  2  cases.  This  causes  large  disparity  in  hr  because  it  is  dividing 
by  a  number  that  is  nearly  equal  to  zero.  This  observation  is  also  made  in  Figure  4-17, 

where  the  highest  heat  flux  cells,  ®/n  >  0.8,  have  a  large  discrepancy  in  heat  transfer 

Vmax 


coefficient. 
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Alpha  =  90,  Mach  =  5 


Figure  4-16  -  Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations  in 
wall  boundary  conditions  at  Mach  5  and 
a  =  90°. 


Figure  4-16  shows  a  break  down  in  predictive  capability  for  a  high  angle  of  attack. 
This  could  be  due  to  grid  dependency  or  complexity  of  the  flow-field,  such  as  highly 
turbulent  recirculation  zones  in  front  of  and  behind  the  airfoil. 


Alpha  =  0,  Mach  =  2 


Variation  in  hr  (%) 


Alpha  =  0,  Mach  =  5 


Figure  4-17  -Comparison  of  ®/n  as  a  Figure  4-18  -  Comparison  of  ®/n  as  a 

vmax  vmax 

function  of  %  variation  in  recovery  function  of  %  variation  in  recovery 
temperature  based  heat  transfer  coefficient  temperature  based  heat  transfer  coefficient 
at  Mach  2  and  a  =  0°.  at  Mach  5  and  a  =  0°. 
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The  large  differences  in  heat  transfer  coefficient  at  the  leading  edge  are  not  present 
in  the  Mach  5  analyses,  where  the  wall  temperature  is  not  nearly  as  close  to  the  recovery 

temperature.  In  fact,  at  ®/n  >  0.3,  we  see  a  +15%  and  -20%  spread  in  difference  in  heat 

Vmax 

transfer  coefficient,  illustrated  in  Figure  4-18. 


Alpha  =  90,  Mach  =  5 
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Figure  4-19  -  ^ / n  as  a  function  of  %  Figure  4-20  -  ®/n  as  a  function  of  % 

vmax  Vmax 

variation  in  recovery  temperature  based  heat  variation  in  recovery  temperature  based  heat 
transfer  coefficient  at  Mach  5  and  a  =  45°.  transfer  coefficient  at  Mach  5  and  a  =  90°. 


The  wall  temperature  dependency  of  the  heat  transfer  coefficient  begin  to  look 
worse  as  the  angle  of  attack  increases  to  45°,  where  the  flow-field  becomes  much  more 
complex,  and  will  likely  be  driven  by  transient  effects.  At  a  —  45°,  Figure  4-19,  the  variation 

in  heat  transfer  coefficient  is  good  for  ®/n  =  1,  where  we  observe  a  ~  ±  2%  spread  in 

Vmax 

heat  transfer  coefficient  variation.  As  seen  in  Figure  4-20,  as  a  ->  90°,  the  difference  band 
becomes  atrocious  and  the  heat  transfer  coefficient  cannot  be  trusted. 

As  mentioned  above,  the  flow-field  becomes  very  complex  at  a  >  45°  due  to  large 
amounts  of  separation.  Figure  4-21  shows  the  separation  occurring  behind  the  airfoil  at  a 
45°  angle  of  attack  and  Mach  5  flow.  Figure  4-22  shows  the  massive  amount  of  separation 
behind  the  airfoil  at  90°  angle  of  attack.  Here,  the  wake  shows  an  unphysical  flow-field.  The 
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dual  recirculation  zones  are  an  artifact  of  the  steady  state  analysis;  in  steady  state,  there 
must  be  symmetry  here  to  balance  the  pressure  in  each  eddy.  In  a  real  flow,  this  would  be 
an  alternating  eddy  from  one  side  to  the  other,  which  causes  a  pressure  misbalance. 


E 


Figure  4-21  -  Separated  flow  behind  the 
airfoil  at  Mach  5  and  =  45° . 


2 

x(m) 


Figure  4-22  -  Separated  flow  behind  the 
airfoil  at  Mach  5  and  =  90° . 


The  large  disparities  in  heat  transfer  coefficient  on  the  windward  side  of  the  Mach  2 
analyses  are  the  result  of  the  wall  temperatures  being  close  to  the  recovery  temperature  on 
the  windward  side.  Figures  4-23  and  4-24  show  the  wall  temperatures  as  a  function  of  path 
length  for  two  selected  flight  conditions.  The  recovery  temperatures,  using  an  effective 
Prandtl  number  of  0.77,  are: 

Tr  =  226.65  K  (l  +  V077  -  22) 

Tr  =  392.8  K\m=2 

and 

Tr  =  226.65  K  (l  +  V077— -  *52) 

Tr  =  1265.3  K\m=5 
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Alpha  =  0,  Mach  =  2 


Alpha  =  0,  Mach  =  5 


Figure  4-23  -  Wall  temperature  as  a  function 
of  path  length  for  Mach  2  and  a  =  0°. 


Figure  4-24  -  Wall  temperature  as  a  function 
of  path  length  for  Mach  5  and  a  —  0°. 


The  above  calculations  can  be  compared  to  Figures  4-23  and  4-24  to  see  why  the 
heat  transfer  coefficients  in  the  Mach  2  cases  have  large  differences  at  the  leading  edge, 
while  the  Mach  5  cases  do  not.  In  Figure  4-23,  the  wall  temperature  at  a  path  length  of  0  m 
is  nearly  identical  to  the  calculated  recovery  temperature  from  above;  this  closeness  in  wall 
temperature  gives  the  effect  of  dividing  by  zero,  and  causes  large  disparity  in  the  calculated 
heat  transfer  coefficient.  In  Figure  4-24,  we  see  that  the  wall  temperature  in  the  same 
location  is  nearly  200  K  below  the  recovery  temperature. 

The  reason  for  the  larger  difference  in  wall  temperature  at  the  windward  side  and 
the  recovery  temperature  for  the  REWT  and  conducting  walls  is  that  these  walls  re-radiate 
heat.  From  the  radiation  heat  transfer  equation,  where  s  is  the  emissivity,  a  is  the  Stefan- 
Boltzman  constant, 

q  =  ecjA(T*-T*)  (4-7) 

We  see  that  the  radiation  heat  transfer  increases  with  the  difference  of  the  fourth  power  of 
the  temperatures.  Thus,  the  wall  temperature  drops  more  significantly  the  higher  the 
temperature  disparity  between  the  free  stream  and  the  wall  temperature. 


55 


Alpha  =  90,  Mach  =  5 


Figure  4-25  -  Wall  temperature  as  a  function 
of  path  length  for  Mach  5  and  =  90°  . 

Figure  4-25  shows  how  differently  the  CFD  software  is  predicting  the  wall 
temperatures  for  the  different  wall  boundary  conditions  at  a  —  90°  and  Mach  5.  These 
differences  in  wall  temperature  illustrate  the  complexity  of  the  flow-field;  as  stated 
previously,  this  case  is  probably  driven  by  transient  effects.  Thus,  it  is  unrealistic  to  run  this 
case  in  steady  state. 

Although  large  variations  in  the  heat  transfer  coefficient  with  respect  to  wall 
temperature  are  observed  in  zones  of  recirculation  and  at  high  angles  of  attack,  we  can  use 
this  method  with  a  moderate  level  of  confidence  for  the  trajectory  in  Chapter  5.  For  the  vast 
majority  of  time,  the  vehicle  travels  at  less  than  45°  angle  of  attack.  Even  during  the  rocket- 
back  maneuver,  the  time  spent  at  these  high  angles  of  attack  is  small,  so  the  temperature 
rise  is  minimal  when  a  heat  transfer  coefficient  is  applied  to  a  material  model  for  a  small 
amount  of  time. 

There  are  regions  of  separation,  even  when  the  vehicle  flies  at  zero  angle  of  attack, 
such  as  in  the  rear  of  the  vehicle,  as  seen  in  Figure  4-26,  where  there  are  rocket  engines 
present.  Thus,  the  separation  induced  variation  in  the  heat  transfer  coefficient  is 
insignificant  in  comparison  to  the  plume  effect,  which  is  not  being  modeled  in  this  study.  It 
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is  likely  that  the  material  selection  in  this  area  will  be  sized  by  the  radiative  heating  from 
the  plume.  The  details  of  this  applied  problem  to  a  3-D  geometry  will  be  discussed  in 
Chapter  5. 


Figure  4-26  -  Illustration  of  the  separation  and  recirculation  around  the  rear  of  the  vehicle. 
The  streamlines  are  plotted  on  the  center  plane.  30%  translucency  is  applied  to  the  center 
plane. 


The  information  provided  in  this  chapter  allows  us  to  move  forward  with  the 
transient  thermal  protection  system  analysis,  covered  in  Chapter  5.  Although  there  are 
some  large  differences  in  heat  transfer  coefficient  in  low  heat  flux  regions,  the  heat  transfer 
coefficient  calculated  by  the  free  stream  temperature  is  much  larger.  For  example,  the 
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largest  variation  in  heat  transfer  coefficient  on  the  200  A" wall  from  Study  B  is  108%,  but 
this  is  in  a  low  heat  flux  region:  ®/n  =  0.03.  Conversely,  for  a  heat  transfer  coefficient 

calculated  using  free  stream  temperature  instead  of  recovery  temperature  on  the  same  wall 
boundary  condition,  has  a  variation  of  2,804%  at  a  ®/n  =1.  Clearly,  we  see  a  large 

Vmax 

dependency  on  wall  temperature  when  free  stream  temperature  is  used  instead  of  recovery 
temperature,  which  accounts  for  the  kinetic  energy  of  the  flow. 

Because  the  variations  in  heat  transfer  coefficient  calculated  by  free  stream 
temperature  are  so  high  when  compared  to  the  heat  transfer  coefficient  calculated  by 
recovery  temperature,  we  see  the  need  to  use  the  one  generated  by  the  recovery 
temperature  when  analyzing  the  transient  temperature  profile  of  a  thermal  protection 
system.  This  process  is  carried  out  in  Chapter  5. 
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CHAPTER  5 


THE  APPLIED  PROBLEM 


The  applied  problem  follows  the  process  flow  shown  in  Figure  5-1.  A  trajectory  and 
geometry  are  provided.  These  parameters  are  taken  into  consideration  to  create  a  universal 
grid  that  will  be  used  in  all  points  solved  in  the  trajectory.  This  universal  grid,  along  with 
the  trajectory,  fed  into  the  external  aerothermal  analysis,  which  is  comprised  of  52  steady 
state  analyses  performed  in  Ansys  Fluent.  These  steady  state  analyses  then  are  used  to 
calculate  recovery  temperature  based  heat  transfer  coefficients  in  the  manner  discussed  in 
Chapters  2  and  4.  These  heat  transfer  coefficients,  along  with  a  material  model  of  a  thermal 
protection  system,  are  applied  to  a  TPS  heat  transfer  analysis,  which  is  a  transient 
calculation  performed  in  ABAQUS. 


Geometry 

Universal 

- * 

External 

Grid 

- > 

Aerothermal 

Analysis 

Trajectory 


TPS  Heat 
Transfer 
Analysis 


Figure  5-1  -  Schematic  of  the  solution  process  used  in  the  applied  problem. 
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A  Rocket-back  Trajectory 

The  rocket-back  trajectory  consists  of  a  rocket-back  maneuver,  where  the  vehicle 
flips  its  tail  into  the  wind,  and  uses  the  main  rocket  engine  to  decelerate  the  booster  after 
second-stage  separation.  This  rocket-back  maneuver  is  apparent  in  Figure  5-2,  which 
graphically  shows  the  vehicle's  altitude  and  angle  of  attack  over  the  entire  mission.  Figure 
5-3  shows  the  Mach  number  and  angle  of  attack.  In  this  trajectory,  the  angle  of  attack  peaks 
just  over  150°  and  the  maximum  altitude  is  slightly  higher  than  300,000  ft.  This  trajectory 
is  conceptual  and  derived  for  the  purpose  of  the  present  academic  research.  It  is  possible 
that  this  trajectory  may  not  be  physically  realistic. 


Figure  5-2  -  Angle  of  attack  and  geodetic  Figure  5-3  -  Mach  number  and  dynamic 
altitude  of  the  given  trajectory.  The  rocket-  pressure  of  the  given  trajectory, 

back  maneuver  starts  at  153  seconds  and 
ends  at  193  seconds 


In  order  to  provide  boundary  conditions  to  the  CFD  code,  the  1976  US  Standard 
Atmosphere  was  used  at  each  trajectory  point  used  in  the  CFD  simulation.  More  about  the 
implementation  of  the  1976  US  Standard  Atmosphere  will  be  discussed  in  the  Problem  Set¬ 
up  and  Solution  Procedure  section. 

Because  the  conceptual  trajectory  is  comprised  of  472  time  steps,  the  trajectory 
used  for  the  steady  state  CFD  analyses  is  coarsened  into  only  52  points  make  the  problem 
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Dynamic  Presure  (lbf/ft2) 


tractable.  Particular  attention  is  given  to  the  time-region  of  the  rocket-back  maneuver,  as 
this  is  where  flight  conditions  change  most  rapidly.  Figures  5-4  and  5-5  describe  the 
trajectory  points  used  in  the  CFD  analysis. 


Figure  5-4  -  Angle  of  attack  of  the  coarsened  Figure  5-5  -  Mach  number  and  dynamic 
trajectory  used  in  the  applied  problem.  pressure  of  the  coarsened  trajectory  used  in 

the  applied  problem. 


Geometry 

The  outer  mold  line  of  the  vehicle  is  a  reusable  booster  system  concept  originally 
designed  by  Northrop  Grumman  Corporation  for  the  United  States  Air  Force  under  a 
contract  with  the  University  of  Dayton  Research  Institute.  This  vehicle  is  designed  as  the 
first  stage  of  a  multi-stage  vehicle  and  is  intended  to  return  to  earth.  Figure  5-6  is  a  CAD 
rendering  of  the  outer  mold  line  of  the  vehicle. 
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Dynamic  Presure  (lbf/ft2) 


Figure  5-6  -  Outer  mold  line  of  the  reusable  Figure  5-7  -  Rear  shot  of  the  outer  mold  line 
booster  system  concept.  of  the  reusable  booster  system.  No  engines 


are  modeled. 

The  vehicle  is  not  modeled  with  an  upper  stage  or  with  engines.  Figure  5-7  shows 
the  rear  of  the  vehicle  with  no  engine  nozzles. 


The  Universal  Mesh 

The  Navier-Stokes  equations  are  solved  for  the  volume  of  the  flow-field.  As  such,  for 
proper  capturing  of  the  boundary  layer,  an  inflation  layer  is  utilized.  This  inflation  layer  is 
sized  to  capture  the  boundary  layer  in  the  entire  trajectory  that  the  CFD  calculations  occur. 
Prandtl's  laminar  boundary  layer  equation  for  a  flat  plate  serves  as  an  estimate  for  the 
boundary  layer  thickness  along  the  vehicle's  surface  (Thompson  505): 


5  = 


\ 


px 

pu0 


(5-1) 


where  S  is  the  boundary  layer  thickness,  p  is  the  viscosity,  x  is  the  distance  from  the  leading 
edge,  p  is  the  density,  and  u0  is  the  free  stream  velocity. 

The  boundary  layer  is  found  to  be  extremely  large  at  the  high  altitudes,  due  to  the 
low  density.  Conversely,  during  the  ascent  phase,  the  velocity  is  high  and  the  density  is  high, 
so  the  boundary  layer  thickness  is  fairly  small.  For  instance,  at  an  altitude  of  188,000  ft,  and 
a  Mach  number  of  5.6,  the  boundary  layer  thickness  at  5  ft  from  the  leading  edge  of  the 
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nose,  which  is  sufficient  enough  distance  from  the  spherical  part  of  the  nose  that  this 


theory  can  be  applied,  is  calculated  to  be: 


= 


3.3  x 

tt  ■  s 


S  =  0.019  ft  =  0.23  in 

for  a  laminar  boundary  layer.  A  turbulent  boundary  layer  is  thicker.  The  inflation  layer's 
overall  thickness  is  approximately6  in,  but  it  expands  exponentially  from  a  very  small 
thickness.  A  cross  section  of  the  grid  near  the  wall  at  the  symmetry  plane  is  shown  in  Figure 
5-8  and  the  inflation  layer  is  shown  in  more  detail  in  Figure  5-9.  The  corresponding  vector 
plot  for  the  boundary  layer,  Figure  5-11,  in  the  region  of  the  above  calculation  shows 
sufficient  capture  of  the  boundary  layer;  there  are  about  six  cells  in  the  boundary  layer. 

As  discussed  previously,  the  heating  at  the  leading  edge  is  driven  by  the  amount  of 
kinetic  energy  in  the  flow.  Therefore,  capturing  the  boundary  layer  is  not  as  critical  in  this 
region.  In  addition,  at  the  leading  edge,  the  steepest  velocity  gradient  is  parallel  to  the 
direction  of  flow,  not  normal  to  the  direction  of  flow  as  with  when  the  flow  is  parallel  with 
the  wall.  The  estimated  boundary  layer  thicknesses  at  5  ft  and  100  ft  from  the  leading  edge 
as  a  function  of  time  during  the  rocket-back  maneuver  is  shown  in  Figure  5-10. 


Figure  5-8  -  The  computational  grid  near  the  wall  of  the  reusable  booster  system  concept. 
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To  create  a  grid  applicable  to  many  angles  of  attack,  the  surfaces  denoting  the 
vehicle  are  placed  inside  a  spherical  volume.  The  spherical  region  that  contains  the  wall 
boundary  is  dubbed  the  dynamic  zone.  The  dynamic  zone  is  rotated  based  on  the  angle  of 
attack  and  bank  angle  of  the  supplied  trajectory.  This  sphere  exists  in  a  cut-out  of  the  far- 
field  domain;  this  region  is  named  the  stationary  zone.  The  interface  of  the  dynamic  and 
stationary  zones  is  a  non-conformal  interface,  which  requires  interpolation  to  pass 
information  between  the  two  regions;  the  interpolation  is  a  source  of  some  error  in  the 
solution.  Pictures  of  the  spherical  region  and  the  entire  computational  grid  are  shown  in 
Figures  5-12  and  5-13,  and  a  view  of  the  non-conformal  interface  is  shown  in  Figure  5-14. 


Figure  5-9  -  Detail  of  the  inflation  layer.  Figure  5-10  -  Estimate  of  the  boundary  layer 

thickness  using  Prandtl's  boundary-layer 
equation. 
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Figure  5-11  -  Vector  plot  of  the  boundary 
layer  thickness  at  Mach  5.6,  altitude  of 
188,000  ft  and,  5  ft  distance  from  the 
leading  edge  of  the  nose. 


Figure  5-12  -  Cross  sectional  view  of  the 
computational  domain  showing  the  dynamic 
zone  (green  spherical  region)  and  the 
stationary  zone  (beige  region). 


Figure  5-13  -  Three-dimensional  sectioned  Figure  5-14  -  Non  conformal  interface, 
view  of  the  dynamic  zone  inside  the 
stationary  zone. 


The  mesh  metrics  are  as  follows.  There  are  26,245,437  volume  elements,  6,083,498 
nodes  in  the  volume,  and  469,396  faces  on  the  wall.  The  maximum  cell  skewness  is  89.98%, 
well  within  the  limit  of  Fluent's  double  precision  solver.  The  maximum  aspect  ratio  is  137.6. 
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Problem  Set-up  and  Solution  Procedure 

With  the  exception  of  the  boundary  conditions,  the  problem  is  set  up  identically  to 
the  one  in  Chapter  4.  Because  the  machine  used  to  solve  the  applied  problem  has  a  Linux 
operating  system  using  a  batch  scheduling  system,  there  is  no  graphical  user  interface  or 
interaction  to  the  solver.  Thus,  journal  files  are  used  to  set  up  the  boundary  conditions, 
rotate  the  mesh,  initialize  the  flow-field,  run  the  calculation,  save  the  case  and  data  files,  and 
export  any  needed  data.  These  journal  files  are  created  in  MATLAB,  which  reference  the 
given  trajectory  in  Excel  format.  This  process  does  not  require  setting  up  unique  case  and 
data  files  for  each  boundary  condition. 

Computational  Resources 

The  steady  state  calculations  were  performed  on  AFRL's  Raptor,  a  Cray  XE6  with  2.4 
GHz  core  speed  and  16  cores  per  node.  The  CPUs  are  AMD  Opteron  64-bit  processors  and 
have  32  GB  of  memory  per  node.  It  uses  a  Cray  Gemini  interconnect  type.  Each  case  file, 
which  contains  the  mesh,  boundary  conditions,  and  solver  settings,  is  about  1.7  GB  of  data. 
Each  data  file,  which  contains  the  solution  data  at  each  element,  is  about  3.9  GB.  Clearly, 
there  is  a  large  amount  of  computational  resource  required  to  solve  each  steady  state 
analysis.  During  the  initial  stages  of  the  computation,  there  was  a  scalability  study 
performed,  the  results  of  which  are  presented  in  Figure  5-15  and  5-16. 
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Number  of  Cores 


Figure  5-15  -  Wall  clock  time  and  CPU  time 
as  a  function  of  number  of  cores  in  the 
scalability  study. 


Figure  5-16  -  Wall  clock  ^  and  CPU  ^  as 

a  function  of  number  of  cores  in  the 
scalability  study. 


The  scalability  study  resulted  in  the  choice  of  200  cores  as  the  optimum  number  of 
processors.  Although  the  200  core  choice  is  not  as  computationally  efficient  as  the  lower 
number  of  cores,  it  is  completed  in  approximately  20  hours.  With  file  transfer  time,  this 
results  in  about  24  hours  per  job.  This  is  a  convenient  time,  as  it  was  fast  enough  that  down 
time  with  the  system  did  not  become  much  of  an  issue,  yet  still  computationally  efficient. 

In  total,  the  total  amount  of  CPU  time  for  the  steady  state  analyses  is  230,924  hours, 
or  26  years,  and  the  total  amount  of  wall  clock  time  is  34  days. 


Energy  Balance 

Before  setting  up  the  transient  TPS  model,  it  must  be  ensured  that  the  heat  transfer 
coefficients  obey  conservation  of  energy.  The  field  of  recovery  temperature  based  heat 
transfer  coefficients  for  one  flight  condition  is  chosen.  This  field  is  then  applied  to  a  surface 
grid  of  the  wall  boundary  only.  This  should  yield  the  same  temperature  profile  as  the  CFD 
results  for  radiation  equilibrium  wall  temperature  at  the  same  flight  conditions.  The  energy 
equation  is  balanced  in  the  following  manner: 
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oe(T*  -  T4)  =  h(Tr  -  Tw)  (5-2) 

where  o  is  the  Stefan-Boltzmann  constant,  a  is  the  emissivity,  Tw  is  the  wall  temperature,  Ts 
is  the  radiative  sink  temperature,  and  Tr  is  the  recovery  temperature.  Indeed,  the  energy 
equation  is  balanced.  Figure  5-17  shows  the  results  from  the  CFD  analysis  using  a  radiation 
equilibrium  wall  temperature  boundary  condition,  and  Figure  5-18  is  the  result  of  the  field 
of  recovery  temperature  based  heat  transfer  coefficients.  These  pictures  are  identical,  which 
proves  the  energy  equation  is  balanced. 


Figure  5-17  -  CFD  results  from  the  radiation 
equilibrium  wall  temperature  boundary 
condition. 


the  balance  of  the  energy  equation,  Equation 
(5-2). 


FEATPS  Model 

Once  the  steady  state  CFD  results  are  complete,  the  heat  transfer  coefficient  is 
calculated  by  exporting  the  heat  flux,  wall  temperature,  and  effective  Prandtl  number  to  a 
tab  delimited  Tecplot  format,  which  is  easily  read  into  MATLAB.  The  data  is  then  exported 
to  a  comma-delimited  format  which  is  easily  read  into  the  Finite  Element  Analysis  (FEA) 
software. 

The  FEA  model  of  the  thermal  protection  system  was  performed  in  ABAQUS.  This 
software  is  a  very  comprehensive  FEA  modeler.  This  software  is  chosen  for  its  ability  to 
solve  a  conduction  heat  transfer  problem  in  three  dimensions  with  a  shell-stack  model  and 
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a  convective  boundary  condition.  This  shell  stack  uses  shell  elements,  which  are  essentially 


surface  elements  of  the  geometry.  Then,  it  applies  integration  points  to  the  elements,  and 
allows  for  composite  layering,  thus  creating  a  shell  stack. 

The  grid  is  extremely  coarse,  as  there  are  not  regions  of  high  gradient  on  the  same 
order  of  magnitude  as  a  CFD  problem,  such  as  in  the  boundary  layer,  or  a  stress  problem, 
such  as  the  area  of  a  stress  riser.  Figure  5-19  shows  the  FEA  grid,  which  can  be  compared  to 
Figure  5-20,  which  is  the  wall  boundary  of  the  CFD  grid.  Note  that  the  CFD  grid  is  especially 
fine  around  the  leading  edges  and  wing  tips. 


z 


Figure  5-19  -  FEA  grid  used  in  the  transient 
heat  transfer  analysis,  which  has  11,912 
shell  elements. 


Figure  5-20  -  CFD  grid  at  the  wall  boundary, 
which  has  469,396  faces  at  the  wall. 


The  material  model  is  a  metallic  thermal  protection  system  concept.  The  material 
thicknesses  are  summarized  in  Table  5-1,  and  the  Material  Properties  are  summarized  in 
Table  5-2.  The  materials  are  modeled  as  constant  density  and  thermal  conductivity  and 
specific  heat  capacity  as  a  function  of  temperature  in  order  to  more  realistically  model  the 
transient  thermal  temperature  profile. 
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Table  5-1  -  Material  thickness  for  the  thermal  protection  system  model. 


Material 

Thickness 

(in) 

Integration 

Points 

Inconel-617 

0.006 

3 

PM-1000 

0.01 

3 

Inconel-617 

0.0015 

3 

PM-1000 

0.01 

3 

Inconel-617 

0.016 

3 

Q-Fiber 

0.5 

3 

Table  5-2  -  Material  thermal  properties  for  the  thermal  protection  system  model. 


Inconel  617 

PM-1000 

Q-Fiber 

T(K) 

u(W\ 

C  (  M 

T(K) 

k  (  w  \ 

Cp 

(  J  ) 

T(K) 

k  (  w  \ 

c  (  ^  \ 

\m  ■  K ) 

p  \k.q  ■  k) 

vm  ■  KJ 

\k.q  ■  KJ 

Vm  ■  KJ 

p  \kq  ■  KJ 

144 

10.52 

385 

294 

12.01 

439 

144 

0.0298 

209 

366 

16.03 

435 

473 

17.01 

502 

311 

0.0298 

643 

478 

18.88 

464 

773 

23.06 

619 

366 

0.0389 

787 

589 

21.12 

490 

1022 

29.98 

720 

422 

0.0476 

932 

700 

23.39 

519 

1272 

37.04 

820 

478 

0.0562 

1025 

811 

25.84 

548 

1472 

41.94 

891 

589 

0.0760 

1079 

922 

28.64 

573 

700 

0.0971 

1111 

1033 

31.41 

602 

811 

0.1195 

1141 

1144 

34.24 

628 

922 

0.1438 

1166 

1366 

39.68 

682 

1033 

0.1731 

1178 

1089 

0.1887 

1182 

1144 

0.2073 

1187 

1255 

0.2422 

1195 

The  model  of  the  thermal  protection  system  is  initialized  at  300  K  and  the  heat 
transfer  coefficient  is  held  constant  through  the  entire  time-step.  The  results  are  reported  at 
the  end  of  each  time-step. 


Results 


The  results  of  the  FEA  are  consistent  with  expectations.  During  the  ascent  phase,  the 
steady  state  temperature  field  over-predicts  temperatures  when  compared  to  the  transient 
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temperature  field,  as  seen  in  Figure  5-21.  This  figure  shows  the  temperature  profile  for  133 
seconds  into  the  flight,  at  Mach  4,  a  4.87°  angle  of  attack,  and  an  altitude  of  130,000  ft  along 
the  symmetry  plane,  y  =  0  m.  The  temperature  profile  has  the  same  characteristic  as  the 
steady  state  temperature  profile,  but  has  obviously  not  reached  equilibrium.  The  peak 
temperature  lies  in  the  region  of  the  nose,  which  corresponds  to  a  path  length  of  0  —  5  m 
and  55  —  60  m.  In  Figures  5-21  through  5-23,  a  path  length  range  of  0  —  30  m  is  the  upper 
surface,  starting  at  the  nose,  and  a  path  length  range  of  31  —  60  m  is  the  lower  surface, 
ending  at  the  nose. 


Angle  of  Attack  4.87  deg,  Mach  4.00,  Altitude  130,000  ft. 


Figure  5-21  -  Temperature  profile  comparison  of  steady  state  CFD  and  transient  TPS 
analysis  for  M  =  4,  a  =  4.87°,  and  z  =  130,000  ft  at  the  symmetry  plane  y  =  0  ft. 
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At  a  mission  time  of  159.5  seconds,  the  vehicle  is  at  a  82.9°  angle  of  attack.  The 
steady  state  CFD  results  reflect  this  high  angle  of  attack  in  its  temperature  profile,  which 
peaks  near  the  nose,  an  area  of  very  low  cross  sectional  area  and  also  still  the  leading  most 
point.  The  temperature  drops  towards  the  rear  of  the  vehicle  at  a  path  length  of  30  m,  as 
seen  in  Figure  5-22.  The  transient  temperature  profile  is  much  different  here.  It  still 
displays  thermal  effects  from  previous  flight  conditions,  most  evidenced  by  the  high 
temperature  at  the  nose  tip  and  the  low  temperature  everywhere  else.  Because  the  vehicle 
spends  very  little  time  at  these  flight  conditions,  there  is  not  enough  heat  transfer  to 
substantially  raise  the  temperature  of  the  thermal  protection  system. 


Angle  of  Attack  82.9  deg,  Mach  6.32,  Altitude  209,000  ft. 


Figure  5-22  -  Temperature  profile  comparison  of  steady  state  CFD  and  transient  TPS 
analysis  for  M  =  6.32,  a  =  82.1°,  and  z  =  209,000  ft  at  the  symmetry  plane  y  =  0  ft. 
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By  the  time  the  vehicle  has  traveled  for  180.7  seconds,  it  has  reached  an  angle  of 
attack  of  151.3°,  an  altitude  of  285,000  ft,  and  has  slowed  to  a  Mach  number  of  3.28.  This  is 
the  peak  angle  of  attack  during  the  rocket-back  maneuver.  Because  the  vehicle  is  now  going 
relatively  slowly,  especially  at  such  a  high  altitude,  the  steady  state  results  under-predict 
the  temperature  profile  of  the  surface  because  the  thermal  protection  has  acted  like  a 
thermal  battery,  storing  the  energy  from  previous  flight  conditions.  This  result  follows 
intuition. 


Angle  of  Attack  151.3  deg,  Mach  3.28,  Altitude  285,000  ft. 


Figure  5-23  -  Temperature  profile  comparison  of  steady  state  CFD  and  transient  TPS 
analysis  for  M  =  3.28,  a  —  151.3°,  and  z  =  285,000  ft  at  the  symmetry  plane  y  -  0  ft. 

Figures  5-24  through  5-29  show  the  temperature  contours  of  the  cases  examined  in 
Figures  5-21  through  5-23.  Note  that  the  temperature  scales  are  not  the  same  in  each  case; 
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this  difference  in  scaling  is  to  allow  for  the  visibility  of  the  temperature  distribution  in  each 


case. 


Appendix  B  summarizes  the  job  name,  mission  time,  integrated  total  surface  heat 
flux,  dynamic  pressure,  angle  of  attack,  and  Mach  number  for  the  CFD  cases  that  were  run, 
which  provided  the  heat  transfer  coefficients  that  are  boundary  conditions  for  the  transient 
model. 


contour  for  M  =  4,  a  =  4.87°,  and 


Figure  5-26  -  Steady  state  temperature 
contour  for  M  =  6.32,  a  =  82.1°,  and 
z  =  209,000  ft. 


Temp(K):  280  300  320  340  360  380  400  420  440  460  480  500  520  540  560 


"igure  5-25  -  Transient  temperature 
contour  for  M  =  4,  a  =  4.87°,  and 
z  =  130,000  ft, 


z 


Temp  (K):  280  300  320  340  360  380  400  420  440  460  480  500  520  540  560  580  600 


Figure  5-27  -  Transient  temperature 
contour  for  M  =  6.32,  a  =  82.1°,  and 
z  =  209,000  ft. 
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Figure  5-28  -  Steady  state  temperature 
contour  for  M  =  3.28,  a  =  151.3°,  and 
z  =  285,000  ft. 


Temp  |K)  290  310  330  350  370  390  410  430  450  470  490  510  530  550  570 


figure  5-29  -  Transient  temperature 
contour  for  M  =  3.28,  a  =  151.3°,  and 
z  =  285,000  ft. 


An  important  point  that  needs  to  be  expressed  is  that  the  transient  analysis 
accounts  for  the  residence  time  of  the  vehicle  at  specific  boundary  conditions  as  well  as  the 
characteristic  heat  transfer  rate.  For  example,  the  “thermal  mass"  can  be  witnessed  in 
Figure  5-23,  where  the  transient  temperature  profile  is  much  larger  than  the  steady  state, 
REWT  temperature  profile.  This  example  shows  that  the  process  follows  reason. 

The  transient  temperature  profile  will  be  affected  by  the  material  selection, 
trajectory,  and  geometry.  In  Chapter  4,  there  were  examinations  of  an  Inconel  617  wall,  a  Ti- 
6A1-4V  wall,  a  radiation-adiabatic  wall,  and  three  isothermal  walls.  We  saw  that  the  same 
analysis  can  be  performed  on  walls  that  have  different  temperature  profiles;  these  different 
temperature  profiles  will  result  from  the  different  material  selections,  different  geometries, 
and  different  trajectories.  Thus,  this  process  can  be  used  on  more  than  the  geometry, 
thermal  protection  system  material  choices,  and  trajectory  used  in  Chapter  5. 
Improvements  to  this  method  and  an  expansion  on  the  significance  of  the  results  are 
presented  in  Chapter  6. 
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CHAPTER  6 


CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 

The  results  presented  in  Chapter  4  allowed  us  to  proceed  with  an  applied  problem, 
which  is  presented  in  Chapter  5.  A  transient  temperature  profile  is  generated  in  Chapter  5; 
this  transient  temperature  profile  can  be  used  to  size  the  thermal  protection  system.  Sizing 
the  thermal  protection  system  means  that  the  materials  can  be  selected,  and  thicknesses 
scaled.  The  process  of  sizing  the  TPS  is  an  iterative  one;  we  would  take  the  generated  heat 
transfer  coefficients  from  the  steady  state  analyses  and  run  numerous  transient  TPS 
analyses  until  we  have  settled  on  a  material  and  thickness  that  satisfies  our  design  criteria. 

The  numerous  heat  transfer  coefficients  from  Chapter  5  can  also  be  used  in  a 
transient  temperature  profile  of  a  structure  coupled  with  a  thermal  protection  system.  The 
heat  transfer  coefficients  serve  as  an  outer  mold  line  boundary  condition,  and  the  thermal 
response  of  the  system  can  be  analyzed,  given  the  heat  generation  rates  of  the  internal 
components  that  generate  heat,  such  as  avionics,  actuators,  batteries,  and  computers. 

Clearly,  we  see  that  the  results  from  this  thesis  mean  that  there  is  a  large 
applicability  of  the  process  used  here,  and  that  the  transient  temperature  profile  can  be 
generated  and  used  instead  of  using  steady  state  results.  Again,  in  Chapter  5,  we  saw  that 
the  steady  state  CFD  results  have  a  large  over-prediction  in  peak  temperature.  For  example, 
in  Figure  5-21,  there  is  approximately  a  250  A" over-prediction  of  wall  temperature  at  the 
nose  of  the  vehicle,  a  path  length  of  2  m.  This  large  over-prediction  is  the  difference 
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between  a  more  expensive  material  and  a  cheaper  one,  or  perhaps  using  a  thinner  TPS, 
which  results  in  weight  and  cost  savings. 

Of  more  importance  is  the  speed  at  which  this  analysis  can  be  carried  out.  In  today's 
rapid  design  environment,  the  use  of  high  fidelity  tools  such  as  CFD  and  FEA  needs  to  take 
place  earlier  in  the  design  phase,  at  the  conceptual  level.  The  use  of  steady  state  analyses  to 
create  a  series  of  heat  transfer  coefficient  fields  coupled  with  a  universal  mesh  allows  this 
rapid  design  to  take  place.  Here,  we  have  been  able  to  demonstrate  that  a  first  pass  design 
at  a  thermal  protection  system  can  be  had  in  slightly  over  a  month's  time,  as  seen  in  Chapter 
5.  This  allows  the  thermal  problems  to  be  attacked  early  in  the  design  phase  instead  of 
pushing  the  problem  back  to  the  final  design  phases,  a  typical  practice  in  current  high  speed 
vehicle  design.  Thus,  the  thermal  management  can  be  integrated  into  the  entire  design  at 
the  conceptual  level.  The  result  is  a  more  accurate  understanding  of  weight,  propellant  mass 
fractions,  cost,  and  material  selection  at  the  conceptual  level. 

Further  Research 

Although  the  recovery  temperature  based  heat  transfer  coefficient  has  a  large  wall 
temperature  dependency  in  low  heat  flux  regions  of  the  wall,  the  method  presented  is 
extremely  useful  to  the  TPS  engineer.  The  method  is  reasonably  quick  to  execute  and  allows 
a  much  greater  level  of  analysis  fidelity  than  previously  available  during  the  early  stages  of 
a  hypersonic  vehicle  design  program.  The  recovery-temperature  based  heat  transfer 
coefficient  more  weakly  dependent  on  wall  temperature  than  free  stream  temperature 
based  heat  transfer  coefficient.  Owing  to  this  fact,  the  use  of  steady  state  CFD  simulations 
can  be  used  in  place  of  the  much  more  time  consuming  transient  CFD  with  conjugate  heat 
transfer  analyses.  The  free  stream  temperature  based  heat  transfer  coefficient  can  be  up  to 


77 


300  times  different  depending  on  the  thermal  boundary  condition  used  at  the  wall.  Clearly, 
accounting  for  the  kinetic  and  thermal  energies  in  the  flow  is  necessary. 

In  an  earlier  discussion,  it  was  noted  that  the  isothermal  walls  have  less  variability 
in  heat  transfer  coefficient  than  the  radiation  equilibrium  wall  temperature  walls.  It  is 
hypothesized  that  this  variability  is  due  to  the  presence  of  a  thermal  gradient  along  the  wall. 
This  hypothesis  is  based  on  the  fact  that  the  flow-field  is  not  greatly  affected  by  the  wall 
temperature  in  high  speed  flow:  the  flow-field  is  dominated  by  the  momentum  of  the  flow 
rather  than  the  temperature  of  the  wall,  which  does  not  hold  for  lower  speed  cases.  The 
energy  is  being  picked  up  near  the  stagnation  point  then  deposited  down-stream;  this 
energy  transfer  is  not  being  accounted  for  in  the  heat  transfer  coefficient  when  thermal 
gradients  are  not  present. 

One  possibility  of  improving  this  predictive  capability  is  to  use  a  local  recovery 
temperature;  that  is,  use  the  velocity  and  temperature  at  the  edge  of  the  boundary  layer. 
This  task  is  much  more  difficult  because  an  algorithm  must  be  developed  to  extract  edge  of 
boundary  layer  data.  Additionally,  this  algorithm  must  handle  zones  with  high  recirculation, 
such  as  the  one  depicted  in  Figure  6-1,  which  is  the  leeward  side  of  the  wing  at  a  90°  angle 
of  attack.  Using  CFD  generated  streamlines  in  two-dimensions  or  path  lines  in  three- 
dimensions  could  serve  as  the  method  for  extraction  of  the  local  recovery  temperature. 

The  use  of  CFD  generated  streamlines  at  the  boundary  layer's  edge  would  allow  for 
extraction  of  pertinent  data,  such  as  the  velocity  and  temperature  of  the  flow  along  the  line 
of  interest.  Additionally,  the  angle  the  streamline  takes  with  the  solid  boundary  may  also  be 
calculated,  which  could  be  used  to  calculate  an  effective  velocity,  or  velocity  projected  to  the 
coordinates  of  the  wall.  This  method  could  more  accurately  predict  the  amount  of  kinetic 
and  thermal  energy  in  the  flow-field  near  that  specific  point  at  the  wall.  Another  approach 
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that  might  be  available  in  the  foreseeable  future,  although  not  cost  effective  today,  is  to  use 


a  Navier-Stokes  flow  solution  coupled  to  a  FEA  TPS  solution. 
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Figure  6-1  -  Turbulent  recirculation  zone  behind  an  airfoil  at  a  90°  angle  of  attack. 


Another  flow  regime  that  the  recovery  temperature  based  heat  transfer  coefficient 
could  be  expanded  to  is  chemically  reacting  hypersonic  flows.  The  recovery  temperature  is 
developed  using  the  ideal  gas  assumption;  an  alternate  method  of  formulation  of  the 
recovery  temperature  for  chemical  dissociation  is  given  by  the  following  equation: 


-h=  1  + 


7-1 


(6-1) 
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where  y  refers  to  the  ratio  of  specific  heats  for  the  undissociated  gas  and  9D  is  the 
characteristic  temperature  of  air  (Liepmann  and  Roshko  351).  This  formulation  of  recovery 
temperature  holds  for  extremely  high  Mach  numbers  because  it  accounts  for  the  heat 
consumed  by  the  chemical  reaction. 

If  the  above  suggestions  are  considered,  we  can  see  a  larger  applicability  of  the 
method  employed  in  this  study,  which  could,  perhaps,  alleviate  the  need  for  transient  CFD 
in  re-entry  and  high  speed  flows  for  many  years.  This  will  result  in  more  rapid  solutions  to 
thermal  problems,  which  can  result  in  faster  design  turn-around. 
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APPENDIX  A 


HEAT  TRANSFER  COEFFICIENT  VARIABILITY  PLOTS  AND  WALL  TEMPERATURE  PLOTS 


Alpha  =15,  Mach  =  2 


Figure  A-l  -  Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations 
in  wall  boundary  conditions  at  Mach  2  and 
a  =  15°. 


Figure  A-2  -  Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations 
in  wall  boundary  conditions  at  Mach  2  and 
a  =  45°. 
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Alpha  =  90,  Mach  =  2  Alpha  =15,  Mach  =  2 


Figure  A-3  -  Comparison  of  %  variation  in 
recovery  temperature  based  heat  transfer 
coefficient  as  a  function  of  path  length  along 
the  surface  of  the  wing  for  each  9  variations 
in  wall  boundary  conditions  at  Mach  2  and 
a  =  90°. 


Figure  A-4  -  ®/n  as  a  function  of  % 

variation  in  recovery  temperature  based 
heat  transfer  coefficient  at  Mach  2  and 
a  =  15°. 


Alpha  =15,  Mach  =5  Alpha  =  45,  Mach  =  2 


Figure  A-5  -  ®/n  as  a  function  of  % 

vmax 

variation  in  recovery  temperature  based 
heat  transfer  coefficient  at  Mach  5  and 
a  =  15°. 


Figure  A-6  -  ®/n  as  a  function  of  % 

Vmax 

variation  in  recovery  temperature  based 
heat  transfer  coefficient  at  Mach  2  and 
a  =  45°. 
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Alpha  =90,  Mach  =  2 


Figure  A- 7  -  ®/n  as  a  function  of  % 

variation  in  recovery  temperature  based 
heat  transfer  coefficient  at  Mach  2  and 
a  =  90°. 


Alpha  =15,  Mach  =  2 


Figure  A-8  -  Wall  temperature  as  a  function 
of  path  length  for  Mach  2  and  =  15° . 


Alpha  =  15,  Mach  =  5 


— —  REWT,  Eps  =  0.9 

- Inconel  617,  Eps  =  0.8 

-•—Inconel  617,  Eps  =  0.9 
Ti-6-4,  Eps  =  0.8 
—  Ti-6-4,  Eps  =  0.9 


Alpha  =  45,  Mach  =  2 


Figure  A-9  -  Wall  temperature  as  a  function  Figure  A-10  -  Wall  temperature  as  a  function 

of  path  length  for  Mach  5  and  =  15°  .  of  path  length  for  Mach  2  and  =  45° . 


85 


Alpha  =  90,  Mach  =  2 
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Figure  A-ll  -  Wall  temperature  as  a  function  Figure  A-12-  Wall  temperature  as  a  function 
of  path  length  for  Mach  5  and  =  45°  .  of  path  length  for  Mach  2  and  =  90° . 
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APPENDIX  B 


WEAK  DEPENDENCE  QUANTIFICATION  STUDY  JOB  CODE  AND  BOUNDARY  CONDITION 

LIST 


Table  B-l  -  Job  code  for  the  weak  dependence  assumption  quantification  study  and  the 
corresponding  boundary  conditions. 


BOUNDARY  CONDITIONS 

Job  Code 

Wall 

Angle  of  Attack 

Mach  Number 

Emissivity 

lim-1-1-1-1 

200  K  Isothermal 

AoA  =  0  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-1-1-2-1 

200  K  Isothermal 

AoA  =  0  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-1-2-1-1 

200  K  Isothermal 

AoA  =  15  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-1-2-2-1 

200  K  Isothermal 

AoA  =  15  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-1-3-1-1 

200  K  Isothermal 

AoA  =  45  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-1-3-2-1 

200  K  Isothermal 

AoA  =  45  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-1-4-1-1 

200  K  Isothermal 

AoA  =  90  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-1-4-2-1 

200  K  Isothermal 

AoA  =  90  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-2-1-1-1 

600  K  Isothermal 

AoA  =  0  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-2-1-2-1 

600  K  Isothermal 

AoA  =  0  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-2-2-1-1 

600  K  Isothermal 

AoA  =  15  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-2-2-2-1 

600  K  Isothermal 

AoA  =  15  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-2-3-1-1 

600  K  Isothermal 

AoA  =  45  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-2-3-2-1 

600  K  Isothermal 

AoA  =  45  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-2-4-1-1 

600  K  Isothermal 

AoA  =  90  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-2-4-2-1 

600  K  Isothermal 

AoA  =  90  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-3-1-1-1 

1000  K  Isothermal 

AoA  =  0  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 
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lim-3-1-2-1 

1000  K  Isothermal 

AoA  =  0  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-3-2-1-1 

1000  K  Isothermal 

AoA  =  15  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-3-2-2-1 

1000  K  Isothermal 

AoA  =  15  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-3-3-1-1 

1000  K  Isothermal 

AoA  =  45  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-3-3-2-1 

1000  K  Isothermal 

AoA  =  45  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-3-4-1-1 

1000  K  Isothermal 

AoA  =  90  deg 

Mach  =  2 

No  Emissivity,  Isothermal  Wall 

lim-3-4-2-1 

1000  K  Isothermal 

AoA  =  90  deg 

Mach  =  5 

No  Emissivity,  Isothermal  Wall 

lim-4-1-1-1 

REWT 

AoA  =  0  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-4-1-1-2 

REWT 

AoA  =  0  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-4-1-2-1 

REWT 

AoA  =  0  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-4-1-2-2 

REWT 

AoA  =  0  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-4-2-1-1 

REWT 

AoA  =15  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-4-2-1-2 

REWT 

AoA  =15  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-4-2-2-1 

REWT 

AoA  =15  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-4-2-2-2 

REWT 

AoA  =15  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-4-3-1-1 

REWT 

AoA  =  45  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-4-3-1-2 

REWT 

AoA  =  45  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-4-3-2-1 

REWT 

AoA  =  45  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-4-3-2-2 

REWT 

AoA  =  45  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-4-4-1-1 

REWT 

AoA  =  90  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-4-4-1-2 

REWT 

AoA  =  90  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-4-4-2-1 

REWT 

AoA  =  90  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-4-4-2-2 

REWT 

AoA  =  90  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-5-1-1-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  0  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-5-1-1-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  0  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-5-1-2-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  0  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-5-1-2-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  0  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-5-2-1-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  15  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-5-2-1-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  15  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-5-2-2-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  15  deg 

Mach  =  5 

Emissivity  =  0.8 
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lim-5-2-2-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  15  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-5-3-1-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  45  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-5-3-1-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  45  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-5-3-2-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  45  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-5-3-2-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  45  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-5-4-1-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  90  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-5-4-1-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  90  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-5-4-2-1 

Radiation  Adiabatic  Ti-6-4 

AoA  =  90  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-5-4-2-2 

Radiation  Adiabatic  Ti-6-4 

AoA  =  90  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-6-1-1-1 

Radiation  Adiabatic  Inconel-617 

AoA  =  0  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-6-1-1-2 

Radiation  Adiabatic  Inconel-617 

AoA  =  0  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-6-1-2-1 

Radiation  Adiabatic  Inconel-617 

AoA  =  0  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-6-1-2-2 

Radiation  Adiabatic  Inconel-617 

AoA  =  0  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-6-2-1-1 

Radiation  Adiabatic  Inconel-617 

AoA  =15  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-6-2-1-2 

Radiation  Adiabatic  Inconel-617 

AoA  =15  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-6-2-2-1 

Radiation  Adiabatic  Inconel-617 

AoA  =15  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-6-2-2-2 

Radiation  Adiabatic  Inconel-617 

AoA  =15  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-6-3-1-1 

Radiation  Adiabatic  Inconel-617 

AoA  =  45  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-6-3-1-2 

Radiation  Adiabatic  Inconel-617 

AoA  =  45  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-6-3-2-1 

Radiation  Adiabatic  Inconel-617 

AoA  =  45  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-6-3-2-2 

Radiation  Adiabatic  Inconel-617 

AoA  =  45  deg 

Mach  =  5 

Emissivity  =  0.9 

lim-6-4-1-1 

Radiation  Adiabatic  Inconel-617 

AoA  =  90  deg 

Mach  =  2 

Emissivity  =  0.8 

lim-6-4-1-2 

Radiation  Adiabatic  Inconel-617 

AoA  =  90  deg 

Mach  =  2 

Emissivity  =  0.9 

lim-6-4-2-1 

Radiation  Adiabatic  Inconel-617 

AoA  =  90  deg 

Mach  =  5 

Emissivity  =  0.8 

lim-6-4-2-2 

Radiation  Adiabatic  Inconel-617 

AoA  =  90  deg 

Mach  =  5 

Emissivity  =  0.9 
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APPENDIX  C 


SUMMARY  OF  STEADY  STATE  3-D  CFD  RUNS  AND  FLIGHT  CONDITIONS 


Table  C-l  -  Job  name,  mission  time,  integrated  total  surface  heat  flux,  Mach  number, 
dynamic  pressure,  and  angle  of  attack  for  each  CFD  analysis  performed  in  the  applied 
problem. 


Job  Name 

Mission 

Time 

(sec) 

Integral 
Total 
Surface 
Heat  Flux 
(Watts) 

Mach 

Number 

Dynamic 

Pressure 

( **.) 
\m  ■  sL) 

ast9500 

95 

448,477 

2.02 

666.5622 

astll507 

115 

1,322,735 

3.00 

303.7527 

astl3300 

133 

1,680,203 

4.00 

92.46446 

astl4800 

148 

969,928 

5.02 

28.39055 

rbtO 

153.45 

727,450 

5.61 

17.98969 

rbtlOO 

154.45 

692,997 

5.72 

16.32025 

rbtl50 

154.95 

684,755 

5.77 

15.53316 

rbt200 

155.45 

672,964 

5.82 

14.76631 

rbt250 

155.95 

657,104 

5.88 

14.03809 

rbt300 

156.45 

738,241 

5.95 

13.43225 

rbt350 

156.95 

892,313 

6.03 

12.82584 

rbt400 

157.45 

1,462,405 

6.10 

12.19808 

rbt450 

157.95 

1,586,430 

6.17 

11.58524 

rbt500 

158.45 

985,463 

6.23 

10.94391 

rbt600 

159.45 

1,153,802 

6.32 

9.588986 

rbt700 

160.45 

693,953 

6.34 

8.180549 

rbt775 

161.2 

551,604 

6.33 

7.157326 

rbt875 

162.2 

451,393 

6.30 

5.944341 

rbt975 

163.2 

369,867 

6.26 

4.894708 

rbtl075 

164.2 

302,757 

6.20 

3.997932 

rbtll75 

165.2 

247,187 

6.14 

3.239218 

rbt!275 

166.2 

197,916 

6.06 

2.604986 
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rbtl375 

167.2 

160,608 

5.97 

2.079738 

rbtl455 

168 

134,293 

5.88 

1.730871 

rbtl550 

168.95 

111,646 

5.77 

1.383836 

rbtl650 

169.95 

91,349 

5.65 

1.090427 

rbtl750 

170.95 

74,741 

5.51 

0.855574 

rbtl850 

171.95 

60,597 

5.36 

0.669147 

rbtl950 

172.95 

49,001 

5.21 

0.5211 

rbt2050 

173.95 

38,823 

5.01 

0.399251 

rbt2150 

174.95 

29,988 

4.77 

0.301348 

rbt2249 

175.94 

23,330 

4.53 

0.227264 

rbt2325 

176.7 

19,086 

4.34 

0.18234 

rbt2425 

177.7 

14,657 

4.08 

0.135936 

rbt2525 

178.7 

11,276 

3.82 

0.100795 

rbt2625 

179.7 

8,697 

3.55 

0.069149 

rbt2725 

180.7 

6,716 

3.28 

0.059024 

rbt2813 

181.58 

5,334 

3.04 

0.050707 

rbt3525 

188.7 

549 

1.20 

0.00792 

rbt3610 

189.55 

486 

1.05 

0.006006 

rbt3700 

190.45 

429 

0.91 

0.004565 

rbt3775 

191.2 

388 

0.84 

0.003838 

rbt3850 

191.95 

465 

0.81 

0.003561 

rbt3907 

192.52 

596 

0.82 

0.003665 

rbt5007 

203.52 

1,056 

0.77 

0.003043 

rbt6207 

215.52 

1,369 

0.76 

0.003545 

rbt9407 

247.52 

7,179 

1.39 

0.02052 

rbtl5807 

311.52 

1,124,225 

2.81 

47.25952 

rbtl6707 

320.52 

1,319,919 

2.86 

228.2062 

rbtl9207 

345.52 

499,918 

2.07 

1029.841 

rbt24407 

397.52 

75,540 

0.70 

179.0888 

rbt41207 

565.52 

129,367 

0.49 

275.2841 
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